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ABSTRACT 

We conduct a series of high-resolution, fully self-consistent dissipationless A^-body simulations to investigate 
the cumulative effect of substructure impacts onto thin disk galaxies in the context of the ACDM paradigm of 
structure formation. Our simulation campaign is based on a hybrid approach combining cosmological sim- 
ulations and controlled numerical experiments. Substructure mass functions, orbital distributions, internal 
structures, and accretion times are culled directly from cosmological simulations of galaxy-sized cold dark 
matter (CDM) halos. We demonstrate that accretions of massive subhalos onto the central regions of host ha- 
los, where the galactic disk resides, since z ~ 1 should be common occurrences. In contrast, extremely few 
satellites in present-day CDM halos are likely to have a significant impact on the disk structure. This is due 
to the fact that massive subhalos with small orbital pericenters that are most capable of strongly perturbing 
the disk become either tidally disrupted or suffer substantial mass loss prior to z = 0. One host halo merger 
history is subsequently used to seed controlled A^-body experiments of repeated satellite encounters with an 
initially-thin galactic disk. These simulations track the effects of six dark matter substructures, with initial 
masses in the range ~ (0.7-2) x 1O 1O M (~ 20-60% of the disk mass), crossing the disk in the past ~ 8 Gyr. 
We demonstrate that these accretion events produce several distinctive morphological signatures in the stellar 
disk including: long-lived, low-surface brightness, ring-like features in the outskirts; significant flares; bars; 
and faint filamentary structures that (spuriously) resemble tidal streams in configuration space. The final distri- 
bution of disk stars exhibits a complex vertical structure that is well-described by a standard "thin-thick" disk 
decomposition, where the "thick" disk component has emerged primarily as a result of the interaction with the 
most massive subhalo. Though our simulation campaign was not designed to elucidate the nature of specific 
Galactic structures, we compare one of the resulting dynamically cold ring-like features in our simulations to 
the Monoceros ring stellar structure in the Milky Way (MW). The comparison shows quantitative agreement 
in both spatial distribution and kinematics, suggesting that such observed complex stellar components may 
arise naturally as disk stars are excited by encounters with CDM substructure. We conclude that satellite-disk 
interactions of the kind expected in ACDM models can induce morphological features in galactic disks that 
are similar to those being discovered in the MW, M31, and in other nearby and distant disk galaxies. These 
results highlight the significant role of CDM substructure in setting the structure of disk galaxies and driving 
galaxy evolution. Upcoming galactic structure surveys and astrometric satellites may be able to distinguish be- 
tween competing cosmological models by testing whether the detailed structure of galactic disks is as excited 
as predicted by the CDM paradigm. 

Subject headings: cosmology: theory — dark matter — galaxies: formation galaxies: dynamics — galaxies: 
structure — methods: numerical 
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In the cold dark matter (CDM) paradigm of structure 
formation, galaxy-sized dark matter halos f orm via the con- 
tinuo us accretion of smaller systems (e.g., iBlumenthal et alj 
1984). Recently, an explosion of data has disclosed un- 
expected kinematic and structural complexity in the Milky 
Way (MW), its neighbors, and distant galaxies. There 
is a growing body of evidence from local star-count sur- 
veys that the MW has, in fact, experienced numerous 
accretion events. Since the original discovery of the 
Sagittarius dwarf galaxy and i ts ass o ciated tidal stream 
(llbataetalJ Il994j lYannv et all l2000tjlvezic et alJ 120001: 
llbata et alJ l2001dTbl : Majewski et al. 2003), at least 5 addi- 



tional, spatially-c oherent structures have been found in and 
around the MW dNewberg et all 120021; iGilmore et all 120021; 
lYannvetafl 12003b iRocha-Pinto et alJ 120041: iMartin et all 
120041: iMartfnez-Delgado et all l2005t iGrillmairl |2006allF 



Grillmair & Dionatos 
iBelokurov et all 120061: 



2006; 



2006; 

Juric et al . 2008). These numbers 
are in general agreement with predictions of dwarf galaxy 
accretion and disruption in the prevailing ACDM cosmo- 



Rocha-Pinto et al] 
I2008h . 
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logical model dBullock et all 1200 U iBullock & Johnstonl 
12005b iBell et all 120071) and die character of the structures 
themselves are similar to those expected in s imulations of 
satellite galaxy disruption dJohnston et al.|[l995l) . 

Among the more intriguing stellar structures discovere d 
in the MW is the Monoceros ring {Rewberg et al. 2002). 
This coherent, low-metallicity feature spans ~ 170° de- 
grees around the Galaxy and lies nea r the disk plane at a 
Galacto-centric distance of ~ 20 kpc dNewberg et alJ 1 2002 1: 
Yannv et alJ[2003t llbata et al.ll2003l:lRocha-Pinto et aLlbooa 
Conn et all 120051) . It is unclear whether this structure is 
yet anot her example of tid a l debris from an a c creted dwarf 
galaxy dYaiinv et alj 120031: iHelmi et alj 120031: iMartin etalJ 
120041: IPenarrubia et al] l2005b or a stellar extension of the 
disk itself ()Ibata et al J2003UHelmi et al.l2003HMomanv et alj 
2004; Moit inho et al.|[2006t iMomanv et al.ll2006l) . 

In addition to revealing a complex assortment of struc- 
ture in the outer disk and halo, photometric surveys are al- 
lowing a more detailed "tomography" o f the MW disk itself 
dNewberg et all 120061: Uuric et al] 12008). Traditional "thin- 
thick" disk decompositions to these data yield thin disk ex- 
ponential scale heights of ~ 250 pc and c orresponding thick 
disk heights typically ~ 2-4 times lar ger dSiegel et al.ll2002h 
iNewberg et al.l 120061: Uuric et alj |2008). However, small but 
significant differences in published scale height determina- 
tions may reflect a more compli cated disk structure th an uni- 
form scale height models allow dNewberg et alj |2006). 

Ongoing studies in the Andromeda galaxy (M31) sug- 
gest that its histor y is at least as rich a nd complex a s that 
of the M W (e g iFerguson et alJ l200l iKaliraiet all 120061: 
Ferguson 2007; Iba ta et al.l 120071, and references therein). 
At least one recent, high-metallicity acc retion event is ev- 
idenced by the "giant stream" structure dlbata et alj 12001a!) 
and Spitzer MIPS imaging has revealed a nearly circular, 
^10 kpc star-forming ring that may have been triggered by 
a past disk interaction dGordon et alJl2006h . Particularly in- 
triguing is the exte nded ~ 40 kpc d isk-like configuration of 
stars around M31 dlbata et alj|2005l) . This extended "disk" 
component is clumpy, but possesses a relatively low veloc- 
ity dispersion, ~ 30kms _1 . It lags behind the rotation ve- 
locity of the standard M31 disk by ~ 40kms _1 and con- 
tains interm ediate age stars similar to those of th e thick disk 
of the MW dBrown et alJl2006b iFaria et al.ll2007l) . Recently, 
IPenarrubia et alJ (l2006h have investigated models where this 
structure is formed by a co-planar, near-circular accretion 
event. Though observations of stellar halos and thick disks 
in more distant galaxies are extremely challenging, deep 
imaging and star-count investigations have also revealed that 
these faint structures are likely ubiquitous (e.g.JSackett et al] 
1994tlDalcanton & Bernstelnl2 002; Zib etti & Fergusorj2004l: 
Dalcanton et alJ 120051 lYoachim & Dalcantonl 120061: iBarthl 
2()()7;lde Jong et alj|2007al) ~ 

The over-abundance of substructure in galaxy-sized dark 
matter halos has received a lot of attention in the literature, 
particularl y regarding the problem of the "missing" galactic 
satell it es dKlypin et aUll999i: iMoore et al.ll 19991: IStoehr et al j 
2002t IZentner & Bullock! 120031: iKazantzidis et all l2004bb 
Kravtsov et alJ2004t Mayer et ~aIll2007l:rStrigari et all2007al) . 
It is important to emphasize that the satellites that will be 
most crucial for altering disk morphologies will be the few 
most massive systems of the entire population that are dis- 
tinctly outside of the regime relevant to the substructure prob- 
lem. Indeed, it can be shown that in the impulsive heating 
approximation the tidal effects of a subhalo population scale 



as dE/dt oc J n(M su \,)M^ uh dM su \„ where n(M su b) andM su b are 
the nu mber density of subhalos and subhalo mass, respec- 
tively (IWhitd l2000). Since cosmological simulations predict 
that the mass function of CDM halos can be d escribed by a 
power law, n(M^b) cx Af~%, with a « 1.8-1.9 (Ghig na et aD 
119981: iGao et ai1l2004l) . the dynamical effects on stellar disks 
will be dominated by the most massive substructures expected 
to be luminous. Therefore, of concern are not the "missing" 
or "dark" subhalos, but rather the population LMC and SMC- 
sized objects that fell into the disk over the past ~ 8 Gyr and 
that are now (presumably) dispers ed throughout the Galaxy 
as part of its ste llar halo (e.g., Bullock & Johnstonl 120051 : 
lAbadi et al.ll2006h . 

In the present paper we examine the morphological signa- 
tures that halo substructure induces in thin galactic disks in the 
context of the CDM model of structure formation. Rather than 
focus on the stellar material liberated from accreted dwarf 
galaxies, we focus on the effects of satellite halo impacts 
on disk galaxies themselves. Specifically, we bombard an 
initially-thin galaxy-sized stellar disk with dark matter subha- 
los whose mass functions, orbital distributions, internal struc- 
tures, and accretion times have been extracted directly from 
cosmological simulations of galaxy-sized dark matter halos 
and study the ramifications of such encounters for disk struc- 
ture. We note that our technique was not designed to repro- 
duce any specific galactic structures and it is worth bearing 
this in mind when comparing our results with those of studies 
aimed at producing specific features. 

Disk galaxies have been notoriously difficult to form i n 
ACDM cosmological simulations (e.g., Navarro et alJl!995l) . 
This difficulty in conjunction with the prevalence of disks 
in the universe might point to an issue of fundamental cos- 
mological importance. For example, the dark matter may 
not be cold or the power spectrum may not be hierarchi- 
cal on small scales (e.g.. ISo mmer- Larsen&Dolgovl 120011: 
IZentner & Bullock 2003; Strig ari et al.ll2007bl) . However, an- 
other possibility is that current ACDM cosmological simu- 
lations lack the physics and/or resolution necessary to form 
disks. Indeed, more recent, higher resolution studies with 
new astrophysica l inputs have been more suc c essful (e.g., 
I Weil et alJ I1998L iThacker & Couchnianl l2001t lAbadi et alJ 
20031: ISommer-Larsen et alj 120031: IGovernato et al.l 12004 
Brook et alJ 120041: iRobertson etaLl l2004t IGovernato et all 
20071) . Nonetheless, the problem of disk galaxy formation 
within a hierarchical cosmology is far from being resolved, 
and any successes are strongly dependent on poorly under- 
stood physical processes that occur below the numerical res- 
olution of contemporary simulations. 

In addition, despite the continuing increase in dynamic 
range, limited numerical resolution renders current cosmolog- 
ical simulations inadequate to address fully problems that en- 
compass the enormous range of dynamical scales involved in 
studies of satellite-disk interactions. Given this fact and the 
uncertainties regarding ab initio formation of disk galaxies, 
we are motivated to explore the interplay between CDM sub- 
structure and galactic disks using controlled A^-body simula- 
tions. Such investigations can have profound implications for 
our current understanding of galaxy formation and evolution. 

Our work is informed considerably by many past 
numerical investigations into the resilience of g alactic 
disks to infallin g satellites dOuinn" & Goodman! IT986fc 
l Ouinn et alJ fl993t IWalker et alJ Il996t iHuang & Carlberd 
19971: ISellwood et alj 119981: I Velazquez & White! | 1999b 
FontetalJ 120011: I Ardi et alT 120031: iGauthier et al.l 120061: 
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Havashi & Chiball2006HRead et al.l2008tlVillalobos & Helrnll 
and we present a more direct comparison to these 
studies in a companion paper (Kazantzidis et al. 2007, Paper 
II). Successes notwithstanding, a significant fraction of these 
earlier numerical investigations had the drawback of not 
being fully self-consistent, modeling various components of 
the primary disk galaxy and sometimes also the satellites as 
rigid potentials, and focusing on experiments with infalling 
systems on nearly circular orbits which are poor approxima- 
tions of the highly eccentric orbits that typify the trajectories 
of CDM substructure. 

Moreover, most of the numerical work published to date 
considered the encounters of single satellites with larger 
multi-com ponent disk gal axies. Notable excep t ions a re the 
studies of lFont et all d2001l) and iGauthier et all d2006l) . who 
utilized the present-day properties of a large ensemble of 
dark matter subhalos in order to investigate the dynamical 
effects of substructure on stellar disks. However, satellite 
popul ations possess vastly different properties at earlier times 
(e.g.JZentner & BuUockll2003HGao et"al]l2004tlZe"ntner et alJ 
l2005ab and investigations based on the z = substructure do 
not account for previous encounters of since disrupted sys- 
tems with the disk. Elucidating the effects of CDM substruc- 
ture on galactic disks clearly requires a more realistic treat- 
ment of the evolution of the subhalo populations. 

Our dissipationless numerical experiments extend and ex- 
pand upon those of earlier studies in at least three major re- 
spects. We adopt for the first time a model in which the sub- 
halo populations are truly representative of those accreted and 
possibly destroyed in the past, instead of the z = surviving 
substructure present in a CDM halo. In particular, our method 
tracks the host halo accretion events since z = 1, and as dis- 
cussed below, this results in a considerably larger number 
of important satellite-disk interactions than previously con- 
sidered. Second, we employ self-consistent satellite models 
whose masses, internal structures, and orbital parameters are 
directly extracted from high-resolution cosmological A^-body 
simulations of galaxy-sized CDM halos. This element of the 
modeling allow us to model the internal properties and im- 
pact parameters of the infalling systems without the need for 
arbitrary choices. 

Finally, the primary disk galaxy models are motivated by 
the prevailing ACDM paradigm of structure formation and 
are derived from explicit distribution functions. In addition, 
they are flexible enough to permit detailed modeling of ac- 
tual disk galaxies, including the MW and M3 1 . In the present 
study we adopt a model that satisfies a broad range of ob- 
servational constraints available for the MW galaxy. The 
self-consistency of the galaxy models in synergy with the 
adopted high-resolution allows us to construct equilibrium N- 
body models of disk galaxies that are as thin as the old, thin 
stellar disk of the MW with an exponential scale height of 
h, ~ 25 ±0.05 kpc (e.g.jKent et al.ll99 U lDehnen & Binnevl 
U998tlWidrow & Dubinskill2005l) . It is worth noting that the 
great majority of past in vestigations (e.g JVelazquez & White! 
1999; Font et all 1200 lb considered A^-body models of disk 
galaxies that were much thicker compared to the old, thin 
disk of the MW. As we demonstrate in Paper II, owing to their 
larger vertical velocity dispersions, thicker disks are relatively 
less affected by encounters with satellites compared to their 
thin counterparts. 

We model the infalling systems as pure dark matter subha- 
los and do not follow the evolution (or deposition) of the stars 
in the accreted systems. Our focus, in this study, is on the evo- 



lution of the stellar material in an initially-thin disk. Though 
our numerical experiments employ a best-fit model for the 
present-day structure of the MW, it is worth emphasizing that 
the adopted simulation campaign is neither designed to fol- 
low the evolution of nor to draw specific conclusions about 
the MW or any other particular system. Given the complex in- 
terplay of effects (e.g., gas cooling, star formation, chemical 
evolution) relevant to the formation and evolution of galac- 
tic disks, our collisionless simulations aim at investigating the 
most generic qualitative features of the evolution of a disk 
galaxy subject to bombardment by CDM substructure. 

We demonstrate that a thin disk component may survive, 
even strongly perturbed, significant bombardment by halo 
substructure and show that a wealth of morphological signa- 
tures are generated in stellar disks in response to these en- 
counters. These include pronounced flares, bars, and a num- 
ber of (sub-dominant) low-surface brightness ring-like fea- 
tures and filamentary structures that arise in and above the 
disk plane. We also show that interactions with infalling satel- 
lites cause stellar disks to develop a complex vertical mor- 
phology that resembles the commonly adopted "thin-thick" 
disk profiles used in the analysis of disk galaxies. The result- 
ing morphological features resemble the complex faint struc- 
tures that are seen in the MW, M31, and in other nearby and 
distant disk galaxies. Quantitative comparison between the 
properties of ring-like features in our simulations and those of 
the Monoceros ring in the MW suggests that such observed 
complex stellar structures may arise naturally as disk stars 
are excited by halo substructure of the kind predicted by the 
ACDM paradigm of structure formation. We reiterate that 
our simulations were not designed to reproduce such features. 
Rather, our campaign was formulated to mimic satellite infall 
in the prevailing model of structure formation and no param- 
eters were tuned to produce these results. The fact that such 
features are excited in our calculations suggests that they are 
of a very general nature, and thus require little fine tuning. All 
of these findings highlight the substantial role of CDM sub- 
structure in setting the observable structure of disk galaxies 
and driving galaxy morphological evolution. 

The outline of this paper is as follows. In Section|2]we use 
cosmological simulations to study substructure properties and 
accretion statistics in four galaxy-sized CDM halos, select tar- 
get subhalos from one host halo merger history for subsequent 
re-simulation, and describe the controlled numerical experi- 
ments to model the impact of these accretion events on the 
structure of an initially-thin galactic disk. Section [3]contains 
results from this simulation campaign pertaining to the global 
morphological evolution that infalling CDM satellites can in- 
duce in galactic stellar disks. Implications and extensions of 
our findings are discussed in Section [4] Finally, in Section [5] 
we summarize our main results and conclusions. Note that 
in the remainder of the paper we use the terms "satellites," 
"subhalos," and "substructures" interchangeably. 

2. NUMERICAL METHODS 

The main goal of the present study is to investigate the mor- 
phological signatures created by the gravitational interaction 
between a population of dark subhalos and a thin galactic 
disk. We begin this section by describing cosmological nu- 
merical simulations of the formation of four galaxy-sized ha- 
los in the concordance ACDM cosmological model. In § 2.2 
we analyze these simulations to quantify the frequency with 
which massive subhalos approach the very central regions of 
their hosts, where the galactic disk presumably lies. We iden- 
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tify one of these host halos, and in § 2.3 we extract from it 
the properties of individual substructures for subsequent con- 
trolled re-simulation. To be specific, we derive the orbital 
parameters, masses, density profiles, and accretion times of 
selected satellites from the cosmological simulations and use 
these properties to initialize controlled numerical simulations 
of subhalo-disk encounters. This allows us to connect the cos- 
mological evolution of a galaxy-sized CDM halo with the evo- 
lution of a thin, galaxy-sized disk. In § 2.4 we introduce the 
primary galaxy models and the numerical methods we em- 
ploy to construct them. Finally, in § 2.5 we describe in detail 
our simulation techniques as well as the controlled numerical 
experiments of the adopted halo merger history. 

2.1. Hierarchical Cosmological Simulations 

The cosmological halos we analyze come from two dif- 
ferent simulations that were performe d using the Adaptive 
Refinement Tre e (ART) TV-body code dKravtsovet all 1 19971; 
lKravtsovl [T999). Both simulations assume a flat ACDM cos- 
mology with Cl m = 1 - fi A = 0.3, fl b = 0.043, h = 0.7 and 
a 8 =0.9. 

The first simulation followed the evolution of three galaxy- 
sized dark matter halos within a cubic volume of 25 /r'Mpc 
on a side using the mu l tiple m ass resolution technique de- 
scribed in Klvpinetal. (2001). We denote this simulation 
"L25" and the halos as "Gj", "G2", and "G3". The mass 
resolution in the vicinity of each system was m p ~ 1.2 x 
10 6 !i~ 1 Mq, corresponding to 1024 3 particles over the entire 
computational box, and the simulation achieved a maximum 
spatial resolution of 191 /i~'pc. Various properties of the L25 
halos and their subs t ructure have been disc usse d previously in 
iKlvpin et al.l (120011) . iKravtsov et all (120041) . and lZentner et al.1 
(20053). 

The second simulation followed the evolution of a single 
galaxy-sized halo in a cubic volume of 20 /i _1 Mpc on a side 
using the same multiple mass technique. We refer to this sim- 
ulation as "L20" and the halo as "G4". The maximum spatial 
resolution for L20 was 153 /r'pc and the mass per particle 
within the region of halo G4 was m p ~ 6.1 x 10 5 /i"'M , cor- 
responding to 1024 3 particles in the entire volume. The L20 
halo has been discussed pr eviously in lPrada e t al. (20061), and 
iGnedin & Kravtsovl (120061) . 

We define virial masses, M v i r , and radii, r vlI , for each host 
halo at z = according to a fixed overdensity A = 337 with 
respect to the mean density of the universe, p, centered on 
the particle with the highest local density. This implies that 
virial mass and radius are related by M v ; r = 47rApr 3 ir /3. Halos 
Gi through G4 have virial masses M v ; r ~ (1 .5, 1 .1 , 1 .1 , 1 .4) x 
10 12 r'Mg, and virial radii r vir ~ (234, 215 , 216, 230) /r'kpc, 
and so they are of sizes appropriate to harbor typical disk 
galaxies such as the MW or M3 1 . Halos G2 and G3 are neigh- 
bors located approximately 427 /i~'kpc from each other in a 
configuration that resembles that of the Local Group. Halo 
Gi is relatively isolated and is ~ 2 Mpc away from the G2-G3 
pair. All three of these halos accrete only a small fraction of 
their final masses and experience no major mergers at z < 1 
(a look-back time of w 8 Gyr), and therefore may reasonably 
host a disk galaxy. Moreover, their merger histories do not 
appear to be extreme compared to the expected spread in halo 
accretion histories measured in simulations with better statis- 
tics ("e.g.. IWechsler et aT1l2002l) . We have chosen Gi as our 
fiducial case for the satellite-disk interaction experiments de- 
scribed in § 12.51 



In order to construct mass accretion histories of the host ha- 
los and orbital trajectories of the accreted subhalos, we iden- 
tified hosts and substructures at a large number of simula- 
tion outputs up to a redshift of approximately z w 10. We 
used 96 available simulation outputs in this interval for L25 
and 148 for L20. We identified halos using a variant of the 
Bound Density Maxima algor ithm which has be en discussed 
extensively in the literature 7 (Kly pin et al.ll 1999b . Virial radii 
are meaningless for subhalos so their radii (and masses) are 
defined using an outer "truncation" radius, r^^. Here we 
adopt r tmnc as the radius at which the logarithmic slope of the 
spherically-averaged density profile of the subhalo exceeds 
a critical value, din p(r)j 'din r\ nmnc = -0.5. This criterion is 
based on the fact that field halo profiles are never shallower 
than this value, and because this definition empirically yields 
a good approximation to the subhalo tidal radius, r t ;d, where 
the background density of the host halo equals the density of 
particles bound to the subhalo. We note here that the spher- 
ical approximation for substructures is fa irly well justified 
dMoore et al.ll2004tlKazantzidis et al.l2005l) and that the outer 
profiles of subhalos typically fall off with dlnp/dlnr < -3 so 
that the bound mass converges well before the formal trunca- 
tion radius is reached. 

The proce dure for tracking subh alo orbital trajectories was 
developed in Kravts ov" et al.l (12004) and we briefly summarize 
it here for completeness. The progenitor for each subhalo at a 
given timestep is selected using the 25% most bound particles. 
We search a previous output timestep for the halo with the 
highest common fraction of these particles and declare this 
system to be the primary progenitor. Many examples of such 
tracks for the L25 halos and a more c omplete description of 
our methods and their test s are given in Kra vtsov et al.l (2004) 
and lZentner et al.l (l2005bl) . 

2.2. Encounters of CDM Substructures with Galactic Disks 

We begin by demonstrating that the z = properties of sub- 
halo populations are strongly biased relative to the history 
of central encounters since a redshift z = 1 . Figure Q] shows 
mass functions and orbital eccentricity distributions of sub- 
structure populations within the halos G1-G4. For this pre- 
sentation, we have scaled the masses and orbital radii of all 
subhalos (M — > f m M and r — > f r r) in order to normalize them 
to a single host halo mass and radius: f m = M w - U /Mk and 
f r = r vk /R h . Here, R h = 244.5 kpc and M h = 7.35 x 10 n M Q 
are the dark halo tidal radius and the dark halo mass interior 
to this radius of the primary disk galaxy model, which we 
use in the controlled satellite-disk encounter simulations de- 
scribed in §2.5. The mass of the disk in this galaxy model is 
M disk = 3.53 x 1O 1O M w 0.048M,,, and we use this number 
to characterize the importance of disk impacts in the figures 
below. Detailed description of the primary disk galaxy model 
will be given in §2.4. 

The left panel of Figure Q] shows cumulative mass func- 
tions of substructures with different orbital pericenters that 
approach the central regions of their hosts, where the galactic 
disk presumably lies. The mass functions are averages over all 
four galaxy-sized host halos in the cosmological simulations 
that we analyze. Satellite masses are measured at the epoch 
closest to when the subhalo first crossed inward a (scaled) 
infall radius of r^ = 50 kpc from the host halo center. We 
define this to be the epoch that our controlled simulations ini- 

7 For completeness, we used rn R i = 5 /i~'kpc and rfl n d = 2.5 /i~'kpc as the 
search radii in the L25 and L20 simulations, respectively. 
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FIG. 1. — Left: Cumulative mass functions of subhalos that approach the central regions of their hosts averaged over four galaxy-sized halos formed in the 
ACDM cosmology. Upper thin lines present results for objects that have had pericenter crossings within 10 (dotted), 20 (solid), and 40 kpc (dashed) since a 
redshift z = 1. The masses are defined at the simulation output time nearest to the first inward crossing of a sphere of radius 50 kpc (in the rescaled units of the 
controlled simulations) and the pericenters are computed from interpolated orbits between timesteps of the simulation outputs. The lower thick lines correspond 
to the (surviving) subhalos identified at z = with pericenters of r pcl j < 40 kpc (dashed) and r pcr i < 20 kpc (solid). These pericenters are also estimates based on 
the orbit of a test particle in a static NFW potential whose properties match those of the host CDM halo at z = 0. CDM models predict several close encounters 
of massive subhalos with the galactic disk since z = 1 ■ In contrast, extremely few satellites in present-day subhalo populations are simultaneously massive 
enough and characterized by small orbital pericenters to have a profound effect on a galactic disk. Right: The distribution of orbital pericentric-to-apocentric 
distance ratios, R = rp er i/r a p , for subhalos with M su b > 0.05Mdi s k. Line types for thin lines are as in the left panel. The thick dot-dashed line corresponds to the 
eccentricity distribution for all substructures at z = 0. Objects that have penetrated deep enough to encounter a disk are on more radial orbits compared to those 
that belong to t he pr esent-day subhalo population. Note that in both panels, all units are relative to the re-scaled units of the controlled satellite-disk encounter 
simulations of $ 12.51 



tiate. The radius rj n f is used to select infalling satellites that 
are likely to have a substantial effect on the disk structure. 
Subhalos that remain on the outskirts of the halo will not af- 
fect the disk as significantly. We refer to satellite properties at 
the simulation output time nearest to the first inward crossing 
of r;„f several times in what follows. The pericenters associ- 
ated with these subhalos are computed from the the orbit of a 
test particle in a static NFW potential whose properties match 
those of the host CDM halo at the time of r m f. 

Many of these satellites are no longer present at z = due 
to disruption, but were present at some earlier epoch. We 
also note that these mass functions are constructed taking into 
account that a single distinct object may pass the host halo 
center multiple times since z = 1 . This panel serves to illus- 
trate that interactions of massive subhalos with the center of 
the host potential where the galactic disk resides since z = 1 
should be common occurrences in the standard ACDM cos- 
mological model. Indeed, in the adopted cosmology, ~ 5 sys- 
tems with masses M su b > 0.2Mdi s k cross through the central 
region (r per i < 20 kpc) of a galaxy-sized halo over a ~ 8 Gyr 
period. 

On the other hand, by examining the present-day substruc- 
ture population of a galaxy-sized halo one would be led to 
believe that close encounters of massive subhalos with the 
galactic disk are very rare: on average, ~ 1 system more 
massive than O.lMdisk is expected to be on a potentially dam- 
aging orbit with r per j < 40 kpc at z = 0. The difference be- 
tween the properties of subhalo populations that were ac- 



creted in the past and those corresponding to the surviv- 
ing systems is not surprising. Subhalos on highly eccen- 
tric orbits that pass closest to the host halo center are pre- 
cisely those that are most likely to lose a large percentage 
of their mass or become tidally disrupted after pericenter 
crossing (e.gJZentner & Bullockll2003l: iKravtsov et al.ll2004t 
iGao et ai1l2004t IZentner et alj 12005 all As orbiting substruc- 
tures continuously lose mass, the fraction of massive satellites 
that interact with disks and are likely to have a significant im- 
pact on the disk structure declines with redshift so that few 
remain by z = 0. 

The right panel of Figure[T]presents orbital eccentricity dis- 
tributions (R = r p eii/r ap0 ) for subhalos with M su b > 0.05Mdi s k 
averaged over the same cosmological simulations. As in the 
left panel, results are presented for systems that have had dif- 
ferent pericenter crossings since z = 1. We also show the 
eccentricity distribution for all subhalos identified at z = 0. 
This is consistent with previous studies reporting a median 
apocenter-to-pericenter ratio of r ap n /r pRr i = 6/1 for pr esent- 
day substructure populations (e.g., Ghi gna et al.lll998l) . The 
plot illuminates another bias in the properties of different sub- 
halo populations: objects that penetrate deep into the center of 
the host potential are on fairly radial orbits, with r per i/r apo < 
0.2, compared to the overall z = subhalo population. It is 
also worth noting that circular orbits appear to be far from the 
norm for systems that sink into the central parts of halos and 
interact with galactic disks. 

While not demonstrated explicitly in the panel discussed 
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above, we find that the orbits of subhalos infalling since 
Z ~ 1 are consistent with isotropic infall from a radius of 
50 kpc. It is known that infall from large r radii is not isotropic 
dKnebe et al.ll2004i : IZentner"era l. 2005b). However, given the 
limited statistical leverage of our sample and our neglect of 
baryonic components in shaping substructure orbits we are 
unable to make a more general statement regarding the direc- 
tionality of satellite infall onto the disk. 

2.3. Selection Criteria and Properties of Simulated 
Cosmological Satellites 

The primary aim of this work is to investigate the morpho- 
logical signatures that arise from a ACDM-motivated satel- 
lite accretion history. While in the previous section we have 
addressed the characteristics of entire subhalo populations 
within host halos G1-G4, here we describe the properties of 
the particular set of infalling CDM substructures that we use 
as the basis for the controlled subhalo-disk encounter simu- 
lations below. We choose to focus on the satellite accretion 
history in host halo Gi discussed in § 2.1. 

After computing the orbital tracks of all subhalos relative 
to the center of the primary host halo Gi, we record all satel- 
lites that cross within a (scaled) infall radius of r m f = 50 kpc 
from the host center since z = 1 . For our re-simulation cam- 
paign, we focus on subhalos that are a significant fraction of 
the disk mass, but not more massive than the disk itself. Fig- 
ure Q] illustrates that, on average, a galaxy-sized dark matter 
halo will have ~ 1 object more massive than Mdisk, and ~ 5 
objects more massive than 0.2Mdi s k cross within the central 
20 kpc since z = 1 . In the present study we select satellites 
with 0.2Mdi s k < M su b < Mdisk- The choice of fairly massive 
subhalos is motivated by the fact that they dominate the tidal 
heating rate by substructure. Neglecting extremely massive 
satellites with M su b > Mdisk is conservative in this case, be- 
cause we aim to study the more subtle features about a rela- 
tively well-preserved disk galaxy, similar to the MW or M3 1 . 
Such massive systems may well destroy or considerably alter 
the morphology of the disk component. In addition, very mas- 
sive satellites would add substantial number of their stars to 
the disk thereby significantly affecting its morphology and ap- 
pearance. We investigate the general robustness of thin disks 
to encounters with systems more massive than the ones con- 
sidered here in a forthcoming study. 

The final selection criterion is related to the orbital percen- 
ters of satellites. To be specific, from the sample of selected 
subhalos, we choose systems which cross the central region 
of halo Gj with r per ; < 20 kpc. Small pericenters are essen- 
tial in order for satellites to be regarded as potentially effi- 
cient perturbers. The imposed selection criteria result in six 
accreted substructures, which we denote S1-S6, to simulate 
over a«8 Gyr period. 

One of the advantages of the present study is that the 
adopted subhalo models are comparatively free from assump- 
tions associated with their internal structure. All previous nu- 
merical investigations of satellite-disk interactions employed 
satellite models with poorly-motivated density profiles rang- 
ing from dens e 7-models to low-co ncentration King profiles 
and truncated iNavarro et al.l (fl996^ hereafter NFW) models. 
The contrast between l ow- and high-density satel lites was pri- 
marily emphasized by Huang & Carlbergj (1 19971) who argued 
that the latter survive tidal stripping and are accreted onto the 
disk, whereas the former are incorporated into outer halos. 
Our approach is to extract the exact density structures of cos- 
mological subhalos and use them to construct high-resolution, 



self-consistent satellite models for the controlled simulations 
below. 

Density profiles of selected cosmological satellites are com- 
puted at the simulation output time nearest to when the sub- 
halo first passes within = 50 kpc of the host halo cen- 
ter. Subhalos S1-S6 show no sign of a central density 
core and are quite cuspy down to the minimum resolved 
radius, as expected from previous work dKazantzidis et al.l 
2004b). Individual infalling substructures are m odeled using 
a multi-parameter (a , 0,~f) density profile law (Zhaol [l996t 
iKravtsov et al.l IT998). A truncation in the density profile is 
also introduced to mimic the effects of the host halo tidal field. 
This function allows a range of inner, 7, and outer, 0, slopes 
across a fit scale radius, r s , and includes a sharpness parame- 
ter, a, to control the nature of the transition between the two 
regimes. The normalization of the density profile is set by a 
characteristic inner density, p s . From a practical perspective, 
truncating the profiles sharply for r > r^, where r t ;d denotes 
the tidal radius of the sate llite, results in models th at are out 
of equilibrium. Following iKazantzidis et all (l2004al) . we im- 
plement an exponential cutoff in the profile of each subhalo 
which sets in at the tidal radius and ensures that for r > r t ;d 
the density of self-bound satellite material falls off smoothly 
to very low values. This procedure results in some additional 
bound mass beyond r^, the precise amount of which depends 
upon the adopted model parameters. In the models considered 
here, the mass exterior to r^ is typically only about ~ 2% of 
the bound mass for each subhalo. 

While it is cumbersome to present the detailed fits for each 
selected cosmological satellite here, we report that they ex- 
hibit NFW behavior in their central parts being very well ap- 
proximated by a steep inner cusp of p(r)otr~ l . In contrast, the 
NFW profile does not provide an accurate description of the 
internal structure of subhalos at intermediate a nd large radii 
in ag r eement with previous nu merical studies (Ghig na et alj 
1998; Kazantzidis et alJl2004bl) . In all cases, the asymptotic 
outer slopes are found to be much steeper than -3. This is 
attributable to the fact that tidal stripping from the host halo 
removes the outer parts of the orbiting satellites steepening 
their density profile. The best-fit values vary from = -3.7 
to = -4.2 suggestin g that subhalos are better represented by 
the lHernquistl dl990|) pro file in their outer parts in agreement 
with lGhigna et all (1 19981) . 

We f ollow the procedure outlined in IKazantzidis et alJ 
(2004a) to construct self-consistent, iV-body realizations of 
satellite models. This method is based on calculating the ex- 
act phase-space distribution function (DF) for a given density 
profile through an Abel transform. For the purposes of our 
study, we assume that the DF depends only on the binding 
energy per unit mass, E. Numerical models of infalling sub- 
structures are then readily generated by randomly sampling 
the particle positions and velocities from this unique DF. A 
critical reader may note that satellites in the present study are 
modeled under the assumptions of spherical symmetry and an 
isotropic velocity dispersion tensor even though cosmological 
halos exhibit strong d epartures from spheri cal symmetry and 
velocity isotropy (e.g. JCole & Laceyll 19961) . 

While this is generally correct, it has been shown that mass 
loss processes (tidal stripping and shocking) inside a host po- 
tential drive substructures toward shapes that are nearly spher- 
ical in both configuration an d velocity space dMoore et alJ 
l2004t IKazantzidis et alJ 12005). Therefore, the assumptions 
of spherical symmetry and velocity isotropy for subhalos 
adopted here are well-grounded. Finally, the masses, posi- 
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FIG. 2. — Cumulative mass profiles, M(r), for representative cosmological 
satellites SI and S2. Stars show profiles derived directly from the cosmolog- 
ical simulation of host halo Gj . Solid curves present fits to the mass distri- 
bution obtained according to the method described in the text (§2.3) and are 
plotted from the adopted force resolution (2e sat = 300 pc) outward. In order 
to emphasize the quality of the fit, we plot the quantity M(r)/r on the vertical 
axis rather than M(r). Masses and radii are normalized to the bound mass, 
M;,, and tidal radius, r tl ^, of each subhalo, respectively. For clarity, the mass 
profile corresponding to the lower curve is vertically shifted downward by 
0.5 dex. The adopted procedure ensures that the satellite models employed in 
this study are initialized according to the density structure of subhalos found 
in high-resolution cosmological iV-body simulations. 



tions and velocities of each selected satellite are scaled in 
order to correspond to the same fractions of virial quantities 
Mjvij, r lv j r , and V v ; r , respectively, in both the parent CDM halo 
and in the primary galaxy model. 

For each satellite model, we generate an Af-body realiza- 
tion containing /V sat = 10 6 particles and employ a gravitational 
softening length of e sat = 150 pc. This choice enables us to 
to resolve density structures down to ~ 1 % of the tidal radius 
of the simulated systems. The adopted particle numbers and 
force resolution are su fficient to resolve mass l oss processes 
inside a host potential dKazantzidis et al.ll2004bl) . 

In what follows, we examine whether the generous al- 
lowance for density profile shapes we considered here can 
give an accurate representation of the density structure of cos- 
mological subhalos. Figure [2] presents cumulative mass pro- 
files, M(r), of characteristic cosmological satellites S 1 and S2 
along with representations of their internal structure obtained 
using the density laws discussed above. Specifically, we plot 
the quantity M{r) /r, where masses and radii are normalized 
to the bound mass and tidal radius of each satellite, respec- 
tively. This choice scales out the gross dependence of the 
mass profile on radius and shows the acceptability of the fit 
better. Figure [2] demonstrates that the satellite models that 
will be used in the numerical experiments of satellite-disk in- 
teractions are initialized according to the internal structure of 
subhalos found in high-resolution cosmological A^-body sim- 
ulations. 

Table [Uprovides a summary of the infall times, orbital pa- 
rameters, and internal structures of cosmological subhalos S 1 - 



S6 that we study in the numerical experiments of satellite-disk 
interactions. Columns (2)-(10) list parameters extracted di- 
rectly from the cosmological simulation of host halo Gi, and 
used as initial conditions in the controlled numerical exper- 
iments. We remind the reader that the structural properties 
for each subhalo were measured at the simulation output time 
closest to the first inward crossing of a sphere of radius 50 kpc 
(in the rescaled units of the controlled simulations). Several 
interesting findings regarding the properties of the simulated 
satellites can be read from the entries in this table. 

First, the 6 different satellite models cover the mass range 
0.21 < M sub /M disk < 0.57, where M disk = 3.53 x 1O 1O M is 
the mass of the stellar disk used in the controlled sat ellite -disk 
encounter simulations that we describe below in § 12.51 This 
mass range corresponds to 7.4 x 10 9 < M su b/M Q <2x 10 10 . 
As we stated earlier, this mass regime is considerably larger 
than the regime of relevance for the missing galactic satellites 
problem (< 10 9 M Q ). As a comparison, the Large Magellanic 
Cloud (LMC) and M32, the most massive dwarf companions 
of MW and M3 1 , have an estimated total present mass of ~ 
2xl0 10 Mo (e.g.. ISchommer et alj|1992t iMastropietro etaT] 
120051) and 2.1 x 10 9 M© (Mate ol 19981) . respectively. The mass 
of satellites S1-S6 corresponds to the upper limit of the mass 
function of observed satellites in the Local Group. 

Another important characteristic of the infalling subhalos is 
their spatial extent. Specifically, they are very extended, with 
r t id > 20 kpc, so of comparable size to the disk itself. This 
fact suggests that the precise pericenter or angular momentum 
vector of the orbit may not be particularly important in driv- 
ing the disk response. We discuss this point in more detail in 
Pap er II. Note tha t with the exceptions of iFont et alj d2001l) 
an d lGauthier et ail d2006l) all previous studies of satellite-disk 
interactions modeled exclusively the compact, baryonic cores 
of the infalling satellites neglecting the more diffuse dark mat- 
ter component. Thus, it remains uncertain whether these ear- 
lier investigations accurately captured the amount of global 
morphological evolution that accreting subhalos can induce 
in cold stellar disks. 

Table Q] also shows that the majority of infalling satellites 
are on fairly eccentric orbits with apocenter-to-pericenter ra- 
tio of r apo /r per i > 6. Interestingly, satellites S4 and S6 appear 
to be on highly radial orbits with r apo /r per i ~ 20. Moreover, 
as we stated earlier, our satellites are consistent with isotropic 
infall but the sample size is insufficient to provide significant 
limits on anisotropic infall. Despite this fact we note the exis- 
tence of an adequate variety of recorded orbits ranging from 
polar (S1,S2) to prograde (S3,S4) to retrograde (S5,S6). 

Finally, the orbital pericenters of the infalling systems mea- 
sured directly in the satellite-disk encounter simulations are, 
in most cases, smaller by a factor of ~ 2 compared to those 
calculated from the cosmological simulation, although there 
are cases where the former are larger (S4 and S6). Several 
reasons may be responsible for these discrepancies. First, the 
presence of the disk in the controlled simulations provides a 
larger central potential compared to the dark halo alone. Sec- 
ondly, the pericenters measured in the cosmological simula- 
tions are computed from the the orbit of a test particle in a 
static NFW potential because the spacing of the simulation 
outputs makes several estimates of pericenters directly from 
orbital tracks rather poor. However, we argue that the global 
disk morphological signatures we report in this study are not 
sensitive to the magnitude of this difference in pericenters. 
The reason is that, as long as a central encounter takes place, 
the subhalos are so extended that their dynamical effects are 
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Properties of the Satellite Models 
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NOTES. — Columns (2)-(10) refer to the properties of cosmological satellites recorded in the simulation output time closest to first crossing within a sphere 
of radius f-j n f = 50 kpc from the center of host halo G] . Col. (1): Satellite model. Col. (2): Redshift at which satellite properties were computed. Col. (3): 
Look-back time in units of Gyr corresponding to the redshift of Col. (2). Col. (4): Satellite distance from the center of the host halo in kpc. Due to the finite 
spacing of simulation outputs this distance differs from 50 kpc. Col. (5): Bound satellite mass in units of the mass of the disk, M disk = 3.53 X 10 10 M Q , used 
in the controlled satellite-disk encounter simulations. Col. (6): Tidal radius in kpc. Col. (7): Maximum circular velocity in kms -1 . Col. (8): Radius at which 
the maximum circular velocity occurs in kpc. Col. (9): Pericentric radius of the orbit in units of the radial scale length of the disk, Rj = 2.82 kpc, used in the 
controlled simulations. Col. (10): Apocentric radius of the orbit in units of Rj. Col. (11): Angle between the angular momentum of the satellite and the initial 
angular momentum of the disk in degrees. This angle is defined so that 8 ~ 90° represents an orbit that starts off near the pole of the disk, 8 < 90° implies a 
prograde orbit, and 8 > 90° corresponds to a retrograde orbit. Col. (12): Pericentric radius of the orbit in units of R^ calculated directly from the controlled 
simulations. This is to be compared with Col. (9). 
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2.4. Primary Disk Galaxy Models 

With the notable exception of [Gauthier et al. (2006j), all 
previous numerical investigations of satellite-disk interac- 
tions used approximate schemes for building TV-body mod- 
els iof_dhk_g^kdxies. Inour study we employ the method 
of IWidrow & Du binski (2005) to construct numerical realiza- 
tions of multi-component disk galaxies consisting of a thin 
stellar disk, a central bulge, and an extended dark matter halo. 
These galaxy models are derived from explicit DFs for each 
component and thus represent axisymmetric, equilibrium so- 
lutions to the coupled collisionless Boltzmann and Poisson 
equations. Because of that, this technique allows us to set up 
thin, equilibrium disk galaxies without instabilities and tran- 
sient effects tha t can be present in other simpler approximate 
schem es (e.g., iHernq uist 1993). The IWidrow & Du binski 
(2005) galaxy models are specified by a large number of pa- 
rameters that permit detailed modeling of actual galaxies such 
as the MW and M31 by fitting to observational data sets. 

For our experiments, we choose their best-fit model MWb, 
which satisfies a broad range of observational constraints 
available for the MW galaxy with simultaneous fits to the in- 
ner and outer rotation curve, bulge velocity dispersion, ver- 
tical force in the solar neighborhood, and total mass at large 
radii. The stellar disk never dominates the rotation curve of 
model MWb. Direct numerical simulations confirm its stabil- 
ity against bar formation for 10 Gyr. In what follows, we give 
a n overview of the parameter s and setup, but refer the reader 
to Widrow & Dubinski (2005) for a detailed discussion. 

The IWidrow & DubinskJ d2005l) disk galaxy models com- 
prise an exponen tial stellar disk, a Hernquist model bulge 
dHernquisl 1 1 9901) and an NFW dark matter halo. The halo 



71 > (1) 
(r/r h ) (l + r/r h ) 

where p/, is a characteristic inner density, and r b denotes the 
scale radius of the profile. For convenience, the inner den- 
sity pi, is expressed in terms of a characteristic velocity 07,, 
ph = al/ATiGr\, with G being the gravitational constant. The 
NFW density profile is formally infinite in extent with a cu- 
mulative mass that diverges as r — > 00. In order to keep the 
total mass finite, it is necessary to truncate the profile. In 
practice, this is accomplished by introducing an energy cut- 
off which is described by the halo tidal radius parameter, £^ 
(0 < £h < 1)- This parameter controls the tidal radius of the 
halo, R/„ which represents the outer edge of the entire system. 
For model building a natural choice for the value of £h would 
be the one for which the resulting R b becomes roughly equal 
to the cosmologically motivated virial radius, r wu . Finally, ar- 
bitrary amount of rotation controlled by the parameter a./, can 
be added to the halo mo del (cei, = 1 12 im plies no rotation). 
The bulge follows the Hernquist ( 1990) density profile 

PH(r) = — — -3 , (2) 

(r/a b ) (l + r/a b ) 

where p b and a/, are a characteristic inner density and scale 
length of the bulge, respectively. In analogy to the treat- 
ment of the halo, the inner density of the bulge p b is ex- 
pressed in terms of a characteristic velocity dispersion a b , 
pb = <j\l 2nGa 2 h . Similarly, the modeling of the stellar bulge is 
complete once the tidal radius parameter, £b, which controls 
the outer edge of the bulge, R b , and the rotation parameter, a b , 
are specified. Finally, the bulge mass is denoted by M b . 

The stellar disks are assumed to be axisymmetric with 
density profiles p L [ = pd(R,z) and DFs which are adopted 
directly from the self-c onsistent disk galaxy models of 
Kuiiken & Dubinski ( 1995). The surface density profile fol- 
lows an exponential distribution in cylindrical radius R, while 
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the vertical structure is mode led by consta nt-thickness self- 
gravitating isothermal sheets ( SpitzeHHH) 

p d (R,z) oc exp (-J^j sech2 ' @) 

where Rd and Zd denote the radial scale length and vertical 
scale height of the disk, respectively. We emphasize that in 
the model formulation the vertical distribution of disk stars is 
described by a sech 2 function. Note that for z > Zd, the sech 2 
law mimics an exponential profile with scale height h z ~Zd/2. 
The construction of the stellar disk is also characterized by 
a truncation radius, R out , and by the parameter, 5R om , which 
controls the sharpness of the truncation. 

The disk phase-space DF is fully determined once the ve- 
locity ellipsoid of the disk is specified. The radial veloc- 
ity dispersion, (Jr(R), is assumed to be exponential with 
<j^(R) = aj, exp(-R/Rd), where ctro denotes the central ra- 
dial velocity dispersion. The disk azimuthal dispersion, 
(JrhiR), is related to ctr(R ) via the epicycle approximation 
dBinnev & Tremaind 19871) . while the vertical velocity disper- 
sion, a z , is set by the requirement that the adopted value of the 
scale height Zd is maintained in the total potential of the galaxy 
model. In all, the density and phase-space DFs for the three 
galactic components are characterized by 14 free parameters 
which may be tuned to fit a wide range of observational data- 
Given a choice of these parameters, the Wi drow & Du binski 
(2005) technique uses an iterative scheme to calculate the 
composite DF and produce equilibrium A^-body models of 
disk galaxies. 

The simulations reported here use Nd = 10 6 particles in 
the disk, N b = 5 x 10 5 in the bulge, and N h = 2 x 10 6 in the 
dark matter halo, and employ a gravitational softening of 
€d = 50 pc, eh = 50 pc, and ef, = 100 pc, respectively. Mass 
and force resolution are adequate to resolve the vertical struc- 
ture of a thin stellar disk and minimize artificial disk heating 
through interactions with the massive halo particles. 

Physical and numerical parameters for each of the compo- 
nents of the adopted primary galaxy model MWb are listed 
in Table 2. The best-fit value of the solar radius and the 
total circular velocity of the galaxy model there are given 
by Rq ~ 8 kpc and V c (Rq) ks 234 km s" 1 , respectively. The 
initial disk scale height is equal to Zd = 400 pc, consistent 
with the observed value for the old, thin stellar disk of the 
MW obtained with a variety o f methods ("e.g.. [Kent et all 
199U iDehnen & Binnevl 119981: iMendez & Guzman! 119981 
Larsen & Humphreys 20031). Obse rvational evidence (e.g., 
Ouillen & GarnettfeOOOtlNordstrom et al.l2004l) indicates that 
the scale height of the thin disk of the MW has not changed 
significantly over the period modeled in this study (z < 1). 
Therefore, it is reasonable to initialize the stellar disk with 
such thickness. 

The Toomre stability parameter of the disk is equal to 
Q = 2.2 at R = 2.5Rd, indicating that the model is stable against 
local non-axisymmetric instabilities. This suggests that any 
significant bar growth during the satellite-disk encounter sim- 
ulations can be assumed to be tidally induced by the infalling 
subhalos. It is also worth emphasizing that the adopted set 
of parameters gives R/, = 244.5 kpc, a value in agreement 
with the nominal virial radius of the MW galaxy as predicted 
by CDM models of cosmological structure formation (e.g., 
iKlvpin et alJl2002h . 

2.5. Description of Satellite-Disk Encounter Simulations 



TABLE 2 

Primary Galaxy Model Parameters 



Component Parameter Value 



Disk: 

M disk 3.53xl0 10 M Q 



Bulge: 



Halo: 





2.82 kpc 




/inn • 
4UU pc 


r> 

«out 


jU KpC 




1 kpc 


CM) 


124.4 kms 


Q(2.5R d ) 


2.2 


Rq 


8 kpc 


N d 


10 6 


e d 


50 pc 


£b 


0.21 




435.7 kms 


"1, 


0.88 kpc 


a b 


0.5 


M b 


1.18 x 10 10 I 


Rb 


3.05 kpc 


N b 


5 x 10 5 


<% 


50 pc 


£ b 


0.1 


&h 


344.7 kms 


rh 


8.82 kpc 


ah 


0.5 



M h 7.35xlO n M 

R,, 244.5 kpc 

N h 2 x 10 6 

e h 100 pc 



All satellite-disk encounter simulations were carried out us- 
ing PKDGRA V, a multi-stepping, parallel, tree A^-body code 
(IStadelll200it) . PKDGRAV uses a spline softening length, 
such that the force is completely Keplerian at twice the quoted 
softening lengths, and multi-stepping based on the local ac- 
celeration of particles. We used an adaptive, kick-drift-kick 
leapfrog integrator with individual particle time steps Af, cho- 
sen according to Af,- < ^(ei/a,-) 1 ' 2 , where e; is the gravita- 
tional softening length of the particle, a, is the value of the 
local acceleration, and rj is a parameter that specifies the size 
of the individual timesteps and consequently the time accu- 
racy of the integration. 

For all simulations, we set the base-timestep to be equal to 
1% of the orbital time of the galaxy model at the disk half- 
mass radius and allowed individual particle timesteps to be at 
most a factor of 2 3l) smaller. In practice, particle timesteps 
were found to be at most a factor of 2 12 smaller than the base 
time-step in the specific application. The number of timesteps 
was sufficient to accurately follow particle orbits down to the 
adopted force resolution of the simulations and resolve the 
vertical frequency of disk stars. The time integration was per- 
formed with high enough accuracy such that the total energy 
was conserved to better than 0. 1 % in all cases. 
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FIG. 3. — Surface brightness maps of disk stars in the simulated accretion history of host halo Gi . The edge-on (upper panels) and face-on (bottom panels) views 
of the disk are displayed in each frame and the color bar in each upper panel indicates the surface brightness limits used to generate the maps. In constructing 
these images, a stellar mass-to-light ratio equal to M*/L = 3 is assumed. Bottom images are 60 kpc on a side, while top images measure 18 kpc by 60 kpc. The 
left panel shows the initial disk assuming that the sequence of satellite-disk interactions initiates at z = 1 ■ The middle and right panels depict the disk after the 
last satellite passage, evolved in isolation for additional ~ 4 Gyr, so that the evolution of disk stars is followed from z = ltoz = 0. In the left and middle panels, 
images are shown to a limit of [i = 26 mag arcsec~ 2 , while the right panel corresponds to a "deeper" surface brightness threshold of = 30 mag arcsec~ 2 . Results 
are presented after centering the disk to its center of mass and rotating it to a new coordinate frame defined by the three principal axes of the total disk inertia 
tensor. Considerable flaring and a wealth of features that they might falsely be identified as tidal streams can be seen in the perturbed disk down to 26 - 30 mag 
arcsec -2 . The existence of non-axisymmetric structures such as extended outer rings and bars after a significant amount of time subsequent to the last accretion 
event confirm their robustness and indicate that axisymmetry in the disk has been destroyed and is not restored at late times. 



We modeled satellite impacts S1-S6 as a sequence of en- 
counters. Starting with subhalo SI we included subsequent 
systems at the epoch when they were recorded in the cosmo- 
logical simulation (Table [TJ. Due to the fact that the simu- 
lated infalling substructures were fairly massive and of com- 
parable size to the disk itself they were not introduced in the 
simulations directly, but rather were grown adiabatically in 
their orbits. This procedure ensures that the disk does not 
suffer substantial perturbations by the sudden pr esence of the 
satell ite and change in potential at its vicinity ( Walk eret al.1 
1996). The mass of each satellite was increased linearly to its 
final value during a timescale that ranges between ~ 150 and 
400 Myr depending on the subhalo mass. The disk structure 
was found to evolve only slightly during the growth period of 
each satellite. 

Accreting subhalos were removed from the controlled sim- 
ulations once they reached their maximum distances from the 
disk after crossing, and thus any given satellite was not per- 
mitted to complete a second orbit. Uniformly, these distances 
were only slightly smaller compared to the starting radii of 
the orbits. The accretion times of simulated subhalos S1-S6 
(Table [U were such that the disk interacted simultaneously 
with the first two satellites SI and S2. We conducted an ad- 
ditional experiment in which subhalo S2 was introduced in 
the controlled simulation after the removal of satellite S 1 and 
confirmed that the disk exhibited the same global evolution 
and morphological signatures as in the standard case. 

In all of the experiments we performed, the satellites lost 
> 80% of their mass after the completion of the first or- 
bit. This justifies our decision to ignore the computation- 
ally costly, but dynamically insignificant, second (subsequent) 
crossing event. Note that only the self-bound core of each 



satellite was extracted from the simulations and not the un- 
bound material, a significant fraction of which remained in 
the region of the disk. Removing the latter would result in 
potential fluctuations that could throw the disk out of equilib- 
rium, and thus interfere with the interpretation of the results. 

In order to save computational time, when time intervals be- 
tween subhalo passages were larger than the timescale needed 
for the disk to relax after the previous interaction, we in- 
troduced the next satellite in the simulation immediately af- 
ter the disk had settled from the previous encounter. Due 
to the complexity of the interaction, we determine this "set- 
tling" timescale empirically by monitoring basic properties 
of the disk structure (e.g., surface density, velocity disper- 
sions, thickness) as a function of time. When these quantities 
stop evolving significantly within radii of interest (r < 6Rd), 
the encounter is deemed complete. Changes of the order of 
5- 10% in disk properties were considered acceptable. We 
note that the limiting radius of 6Rd in chosen because it ini- 
tially contains ~ 98% of the mass of the disk. Typical settling 
timescales correspond to ~ 100-200 Myr after the removal 
of the bound core of each orbiting subhalo. This eliminates 
the computational expense of simulating the disk during the 
quiet intervals between interactions. 

Finally, we stress that prior to commencing the satellite- 
disk encounter simulations, the primary galaxy model was 
tilted so that the disk angular momentum vector was aligned 
with the angular momentum of the host CDM halo Gi. This 
choice is motivated by results of cosmological simulations 
suggesting that the angular moment a of galaxies and their host 
halos tend to be well aligned (e.g jLibeskind et aUl2007t) . 

3. RESULTS 
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In this section we present results regarding the evolution 
of an initially-thin stellar disk subject to a cosmologically- 
motivated subhalo merger history. Recall that we examine 
the disk response to substructure impacts S1-S6, which are 
designed to mimic a reasonable central accretion history for 
a galaxy-sized halo over the past ~ 8 Gyr. The "final" disk 
discussed in the next sections has experienced the S1-S6 en- 
counters and was evolved in isolation for ~ 4 Gyr after the 
last interaction. While this allows for relaxation after the en- 
counters, it also ensures that all of the resultant morphological 
features are long-lived rather than transient. Consequently, 
our results are relevant to systems that exhibit no obvious, on- 
going encounters and have been allowed to relax in isolation. 
We also stress that we do not include the bulge component 
in any of the analysis presented below. The focus is only on 
the evolution of the disk material. We defer any discussion 
regarding the bulge to future work. 

We compute all properties of the disk and show all visual- 
izations of the disk morphology after centering the disk to its 
center of mass and rotating it to a new coordinate frame de- 
fined by the three principal axes of the total disk inertia tensor. 
The motivation behind performing both actions is twofold. 
As discussed in more detail in Paper II, exchange of angu- 
lar momentum between the infalling satellites and the disk 
tilt the disk plane substantially over the merging history and 
cause the disk center of mass to drift from its initial position 
at the origin of the coordinate frame. In the original coordi- 
nate frame, rotation in a tilted disk would appear as vertical 
motion interfering with the interpretation of the results. 

Finally, we underscore that there is nothing exceptional 
about host halo Gi as far as the properties of the subhalo popu- 
lations are concerned. In all other three galaxy-sized halos we 
analyzed, substructures of similar numbers, masses, internal 
structures, orbital parameters, and accretion times were iden- 
tified. Though our simulation program is designed to mimic 
the activity in the inner halo of Gi, the similarity of subhalo 
populations in all four halos suggests that the results presented 
next should be regarded as fairly general. 

3.1. Global Disk Morphology 

Figure [3] depicts the evolution of the global structure of a 
thin stellar disk that experiences a merging history of the type 
expected in ACDM. The left panel shows the initial config- 
uration of disk stars assuming that the sequence of satellite- 
disk interactions initiates at z = 1, while the right and mid- 
dle panels show the final configuration at z = of the same 
stars at two different surface brightness thresholds. Stellar 
particles are color-coded by the projected surface brightness. 
In producing these images, we have assumed a stellar mass- 
to-light ratio of M*/L = 3 in the V-band, as would be ap- 
propriate for an old stellar population with colors relevant 
for thick disks (Dalcanton & Bernstein 2002). This partic- 
ular choice produces a peak edge-on central surface bright- 
ness of /i ~ 21 mag arcsec -2 , similar to values seen in op- 
tical surveys of nearby, edge -on, undisturbed disk galaxies 
dYoachim & Dalcanton! 120061) . A different choice for M*/L 
would simply result in a dimming or brightening of these im- 
ages and subsequent figures accordingly. 

Several interesting morphological signatures of satellite- 
disk interactions are clear from these images. First, a wealth 
of low-surface brightness features have developed both in the 
plane and above the plane of the final disk as a consequence 
of these disturbances and we discuss these in detail in the next 
sections. Particularly intriguing is the presence of a conspic- 



uous flare and non-axisymmetric structures such as extended 
outer rings and bars. Their existence and robustness indicate 
that the axisymmetry of the disk has been destroyed by the 
encounters with the infalling satellites and is not restored at 
late times. 

Moreover, as is evident in the edge-on images of the final 
disk in Figure [3] a high-surface brightness thin disk compo- 
nent remains after the bombardment by CDM substructure. 
We discuss the implication of this finding in the general con- 
text of disk survivability in a ACDM universe in Section|4] It 
is interesting to ask whether or not the bright in-plane struc- 
ture would be recognized as a thin disk component in a stan- 
dard disk decomposition analysis. We present such an analy- 
sis in Figure [4] This figure shows the initial and final surface 
brightness profiles, /i(z), of the stellar disk as a function of 
distance, z, above the disk plane. For this calculation, the 
stellar distribution is observed edge-on. Symbols correspond 
to different projected radii, R. 

The left panel shows the analytic sech 2 surface brightness 
profile (eq. (3)) for the initial, thin stellar disk with a scale 
height of Zd = 0.4 kpc. By construction, there is excellent 
agreement between the analytic profiles and results obtained 
directly from the particle distribution. As we discussed ear- 
lier, Figure [3] demonstrates that the disk is substantially per- 
turbed by the encounters with the infalling satellites. This 
suggests that describing the later stages of the disk evolution 
requires more than the simple, single function that works well 
for the initial disk. Indeed, the right panel of Figure|4]presents 
a two-component, "thin-thick" disk decomposition for the fi- 
nal disk. The decomposition consists of a low-surface bright- 
ness, "thick" disk component with a sech 2 scale height of 
ZtMck =16 kpc and a "thin" disk component with a sech 2 scale 
height of z t hin = 0.6 kpc which is slightly larger compared to 
that of the initial thin disk. 

This simple decomposition provides a reasonably good de- 
scription of the final disk down to fairly low surface bright- 
nesses, /i < 27 mag arcsec" 2 , and for \R\ < 9 kpc. Interest- 
ingly, the edge-on thick disk component is consistent with 
having a surface brightness profile that is almost independent 
of projected radius for R < 9 kpc. At fainter magnitudes, a 
flare and other diffuse structures become important as dis- 
cussed in § 13.21 Beyond R ~ 9 kpc, the adopted decompo- 
sition fails to provide an adequate representation of the disk 
structure. At these radii the disk is significantly flared and its 
structure becomes strongly irregular. The complexity of the 
disk structure at R > 9 kpc suggests that a more complicated 
functional form may be needed to describe the final disk at 
larger radii. 

We note that this decomposition does not constitute a for- 
mal fit but merely provides an acceptable parametrization for 
the final disk structure. The analysis presented above simply 
serves to illustrate that the final disk cannot be represented by 
a single exponential or sech 2 function. Two components are 
clearly required, one of which is significantly thicker than the 
"thin" one. Deriving a formal fit for the disk vertical structure 
will not yield substantial additional information regarding the 
global morphological evolution of the disk, though it will be 
explored in the future. 

Figure|4]shows a number of interesting results related to the 
surface brightness distributions of the final disk. First, these 
distributions become uniformly broader extending at consid- 
erably larger vertical distances from the disk plane compared 
to those of the initial thin disk. This change is caused by the 
encounters with the infalling subhalos and can be directly as- 
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FIG. 4. — Thin-thick disk decomposition analysis. Surface brightness profiles for disk stars as a function of distance above the plane of the disk, z. The stellar 
distribution is observed edge-on and results are shown for both the initial (left panel) and final (right panel) disk. The final disk has experienced the S1-S6 
encounters and was further evolved in isolation for ~ 4 Gyr after the last interaction to ensure that is has reached a relaxed state. Different symbols correspond 
to various projected radii, \R\ = 0,4,6,9, 12 kpc, averaged between the i\R\ projections about R = 0, and within 1 kpc strips centered on each of the projected 
radii. The solid lines in the left panel present the analytic sech 2 surface brightness profiles of the initial thin disk with Zd = 0.4 kpc and demonstrate a very good 
agreement with results obtained directly from the particle distribution. Solid lines in the right panel show a "thin-thick" disk decomposition for the final disk, 
which consists of a "thick" disk component (dashed line) with a sech" scale height of Zthick = 1.6 kpc and a "thin" disk component (dotted lines) with a slightly 
larger scale height (z t hin = 0.6 kpc) compared to that of the initial thin disk. The surface brightness profile of the "thick" disk is nearly independent of projected 
radius for R < 9 kpc explaining the existence of only one dashed line. At small projected radii (\R\ < 6 kpc), the contribution of the thin component to the light 
profile of the stellar disk is dominant for heights z < 1 kpc. As \R\ increases, the situation is reversed with the thick disk component becoming prevalent. This 
simple decomposition provides a reasonable description of the final disk for p < 27 mag arcsec~ 2 and \R\ < 9 kpc. At fainter magnitudes and beyond R ~ 9 kpc, 
the flare and other diffuse structures become important and the adopted decomposition fails to provide an accurate description of the disk structure. 



sociated with the strong flare seen in the final disk which will 
be discussed in § 13.31 Second, the peak of the surface bright- 
ness profiles which corresponds to the light distribution in the 
disk plane (z = 0) drops systematically at all projected radii 
compared to that of the initial disk. This is attributable to the 
fact that stars are efficiently heated above the disk plane by 
the merging subhalos. As a result of this, the number of stars 
in the plane of the disk decreases with obvious consequences 
for the surface brightness distribution. 

As the projected radius increases, the height above the plane 
of the disk at which the thick disk component dominates the 
light profile becomes progressively smaller. This is because 
the surface brightness distribution of the thin disk declines ex- 
ponentially with R, while both the scale height, and the cen- 
tral (projected) surface brightness of the edge-on thick disk 
vary much more slowly with radius. At small projected radii 
(\R\ < 6 kpc) the thin disk component dominates for z< 1 kpc. 
However, this picture is quickly reversed. For \R \ = 9 kpc (tri- 
angles) the thick disk dominates the light profile of the stellar 
distribution for all heights above the disk plane. We note that 
within 9 kpc, approximately ~ 17% of the final stellar mass is 
contained in the thick disk component. 

3.2. Morphological Signatures of Satellite-Disk Encounters 

The results in the previous section demonstrate that in- 
falling CDM substructures play a significant role in setting 
the global morphology of a stellar disk. The face-on view 



of the final disk in Figure [3] (middle and right panels) re- 
veals that satellite accretion can also excite strong bars which 
drive further evolution in the inner disk regions. The primary 
disk galaxy model was constructed to be stable against the 
formation of a bar. Thus, the observed bar growth should 
be regarded as tidally induced by the infalling satellites. In- 
terestingly, the edge-on view of the same panels shows the 
generation of a characteristic "X" shape in the bright central 
disk, a finding also report ed in previous numerical studies of 
satell ite-disk encounters (Wa lker et al.1 1 19961: iGauthier et al.1 
2006). This noticeable feature is often linked to secular evo- 
lution of galaxies driven by the presence of a bar when it 
buckl es as a result of becomi ng unstable to bending modes 
(e.g.. lCombes & Sand ers 1981), and may be associated with 
the prese nce of peanut-shaped bulges observed in many galax- 
ies (e.g., iLutticke et al1l2000h . 

Moreover, while the bright "X"-shaped component is visi- 
ble at high-surface brightness, many other interaction-driven 
signatures appear as low-surface brightness features. The 
deep view of the final edge-on disk in Figure [3] (right panel) 
shows a number of additional filamentary structures at /i > 26 
mag arcsec -2 , and other complex configurations, that develop 
above the disk plane. These structures bear some resemblance 
to tidal streams, but are in fact disk stars that have been grav- 
itationally excited by the subhalo impacts. The final face-on 
disk in the same panel is also significantly more structured at 
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FIG. 5. — Monoceros ring comparison. Circles are data of M giant stars from the Crane et al. 1 2003) study of the Monoceros stream in the direction of the 
Galactic anticenter. Points in each panel correspond to stars in the outer regions (R > 15 kpc) of the final disk. These stars are color-coded according to their 
galacto-centric radius, with dark blue/black at R = 15 kpc, light-blue at R ~ 18 kpc, green at R ~ 22 kpc, and orange/red at R > 25 kpc. The upper-left panel shows 
a face-on view of the disk and the asterisk marks the adopted position of the Sun at (jt, y, z) = (8, 0, 0) kpc. The upper-right, lower-left, and lower-right panels show 
vertical distance of stars above the disk plane, z, heliocentric distance, and heliocentric radial velocity, respectively, as a function of Galactic longitude in degrees, 
/. Diamonds in the lower-right panel correspond to the median line-of-sight velocities in four bins in longitude / = [240 — 210, 210—180, 180- 150, 150— 120]° 
for a specific stream-feature identified in the simulations. The associated bars reflect the line-of-sight velocity dispersion in each bin. Accretion histories of the 
kind expected in ACDM models produce dynamically cold ring-like features around galactic disks that are quantitatively similar to the Monoceros ring in the 
MW. 



faint surface brightness and large radius compared to the ini- 
tial disk and is quit e reminiscent o f the outer disk structure 
seen around M31 dlbata et alj |2005). The same image also re- 
veals a pronounced ring-like feature at fi ~ 26 mag arcsec" 2 . 
The ring is nearly in-plane and is qualitatively similar to the 
Monoceros stream known to extend around the MW (e.g., 
iNewberg et al.l2002trWn"nv et alj2003l:llbata et alj2003l) . We 
emphasize that this ring is a non-transient feature which sur- 
vives for a considerable amount of time after the satellite ac- 
cretion events. Indeed, though this structure is produced after 
the interaction with satellite S2, it remains apparent at the fi- 
nal simulation output some ~ 4 Gyr after the final impact with 
subhalo S6. 

Though our simulations were neither designed to follow 
the evolution of nor to draw specific conclusions about the 
MW or any other particular system, the reminiscence of the 
resultant rings in our simulations to that of the Monoceros 
ring in the Galaxy is suggestive of the possibility that the 
latter may have been generated via an excitation of an an- 
cient disk's stars. In order to check the general validity of 
this scenario Figure [5]presents a more direct comparison be- 



tween the ring-like structures generated via satellite-disk in- 
teractions and the Monocero s ring feature towards the Galac- 
tic anticenter discovered by (Ne wberg et alJl2002h . The cir- 
cles in t his figure cor r espon d to a kinematic study of M giant 
stars by [Cr ane et al. ( 2003), which followed up a study by 
iRocha-Pinto et al.1 ((2P03) to identify M giants ass ociated with 
the rin g. Our specific choice to compare to the ICrane et ail 
d2003l) sample is motivated by the fact that these data span 
the Monoceros stream uniformly over nearly ~ 100° in the 
sky with precise membership criteria and good velocity deter- 
minations. For clarity, we have focused our analysis on the 
distribution of stars at galacto-centric radii larger than 15 kpc. 
These stars are color-coded according to their radial position 
from the center of the disk. Figure [5] shows vertical distance 
above the disk plane (z), heliocentric distance, and heliocen- 
tric radial velocity of the same stars as a function of Galactic 
longitude, I. 

It is worth emphasizing that the stars in our sim ulated data 
are m ore finely sampled compared to those of ICrane et alJ 
(2003). In order to facilitate the comparison, in the lower- 
right panel, we have more precisely identified a stream-feature 
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FIG. 6. — Effect of individual satellites on the global morphological evolution of the disk. Each panel shows edge-on surface brightness maps of disk stars to a 
surface brightness limit of ^ = 30 mag arcsec and measures 18 kpc by 60 kpc. Labels for individual satellite passages from SI to S5 are indicated in the lower 
left-hand corners of each image and, to aid comparison, the initial disk is presented in the upper left panel. Images are constructed after the disk has relaxed 
from each encounter according to the definition in §2.5. The response to subhalo S6 can be seen in Figure[5] A high-surface brightness thin disk component 
is visible in all images. The first satellite passage SI generates a warp in the outer regions of the disk. The second encounter which involves the most massive 
subhalo S2 causes substantial flaring above the disk plane. The remaining satellites perturb the disk structure further even though their effect is differentially less 
significant compared to the initial encounters. In broad terms, most of the morphological signatures have emerged as a result of the interaction with the most 
massive subhalo. 



along four bins in longitude / = [240-210,210- 180, 180- 
150,150- 120]° and associated cuts in helio-centric dis- 
tance dhdio = [14- 12 kpc, 12.5-10 kpc, 12.5-10 kpc, 12.5- 
10 kpc]. The total mass of this stream-feature is ~ 10 8 M© or 
~ 0.3% of the disk mass, M^. Although the total mass of the 
Monoceros stream is still a matter of debate, the aforentioned 
value is consistent with the mass estimates of lYannv et al.l 
(120031) . 

The four diamonds correspond to the median line-of-sight 
velocities and longitudes in each of these angular bins. The 
bars reflect the line-of-sight velocity dispersion in each bin, 
which is approximately a\ os g 40 kms" 1 across the stream. 
We emphasize that using the ICrane et ail (120031) data across 
the first three of these bins we measured a line-of-sight ve- 
locity dispersion of cxi os ~ 39 kms -1 , which is in excellent 
agreement with the corre sponding values ass ociated with the 
simulated data. Note that ICrane et al.1 (12003) estimate a line- 
of-sight velocity dispersion of <ti os = 20 ±4 km s" 1 for the ring, 
which differs significantly from the aforementioned value of 
<7i 0S — 39 kms" 1 . This is because the former was calculated 
using a subset of 53 of the 58 stars in their sample, ob- 
tained after a 2.5tr (±50 kms" 1 ) rejection threshold was ap- 
plied about a third-order polynomial fit to the velocity trend. 
When we perform a similar rejection threshold within each of 
our four angular bins we find tri os ~ 24.5 kms" 1 for our sim- 
ulated ring feature. Therefore, following the same analysis 
methods, our results are in good agreement with the quoted 
veloci ty dispersion for Monoceros stream from ICrane et al.1 
J2003h . 

The general agreement in the properties of rings (spatial 
distribution and kinematics) generated in our simulations with 



those of the Monoceros ring structure from the ICrane et all 
(2003) data is encouraging. This agreement is particularly 
noteworthy because our simulation program did not aim to 
reproduce such a feature. Our simulation campaign has ex- 
plored only one realization of a galaxy-sized halo accretion 
history and we tuned no aspect of the computations to produce 
this agreement. This is a significant point that bears emphasis. 
One may constrast our model of a disk origin for the Mono- 
ceros stream to previous ones that have attempted to explain 
this structure via the accretion and disrupti on of a satellite 
galax y on a nearly-circular, co-planar orbit (Penarrubia et al. 
2005). Neither our model nor the model of iPenarrubia et al.l 
(2005, see Figures 2 & 4 in this reference) can be falsified 
by the extant data; however, our ring was produced as an 
unforeseen bypr o duct o f our simulation campaign while the 
Penarrubia et al] d2005l) ring was produced as part of a pro- 
gram aimed at tuning a satellite accretion event to yield a 
structure similar to the Monoceros ring. As such, a scenario 
of the kind we propose seems a viable mechanism to produce 
such ring-like structure from a stellar disk. 

Given the existence of at least two qualitatively different 
scenarios for the formation of structures like the Monoceros 
ring, it is crucial to be able to test these formation scenarios 
and distinguish between them. In principle, the metallicity of 
stars in the stream may constrain the competing models. The 
stars in an accreted satellite can be generally expected to have 
different metallicity from that of the surrounding disk stars, 
while in our model the metallicity of the stream should be 
comparable to the metallicity of the thick disk stars surround- 
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ing the stream. Unfortunately, there is no definitive conclu- 
sion regarding the metallicity distribution of the Monoceros 
ring stars as the metallicity estimates span a wide range of 
values. Similar uncertainties exist for the characterization of 
the outer t hick disk. 

Indeed, Yann y et ail (120031) report mean metallicities for 
their sam ple of Monoceros s tream stars of [Fe/H] = — 1.6 ± 
0.3, while ICrane eTail f2003) measure a much higher metal- 
licity of [Fe/H] = -0.4 ± 0.3 for the M giants associated with 
the ring. The mo st recent measurements from the SDSS 
dlvezic et al.ll2008l) indicate a metallicity for the Monoceros 
ring of [Fe/H] « — 1, which is quite similar to their own de- 
terminations for the outer, high-latitude disk, [Fe/H] w -0.9. 
Other determinations of the thick disk metallicity distribution 
have shown a dependence on how one differentiates thick disk 
from halo stars, but the t ypical range is [Fe/H] sa -0.7 to 
-L8 dGilmore et alj|1995t IChiba & Beersll2000t iBrown et al] 
2008). The lower ring metallicities reported by Yan ny et al.l 
( 2003 ) are consi stent with the lower range of thick-disk metal- 
licities quoted (Bro wn et al.l 120081) , but would be problem- 
atic for our model if the outer thick disk is truly as metal 
rich as the canonica l thick-disk value, [Fe/H] rs -0.7 to -1 
dGilmore et al.|[l99l llvezic et al.ll2008l) . Taking this canoni- 
cal thick-disk valu e, th e higher ring metalli cities estimated by 
ICrane et all d2003l) and llvezic etal.1 d2008l) would seem to be 
more in line with our scenario. Overall, the use of metallicity 
measurements as a robust constraint on our model would re- 
quire refining the observational measurements for the Mono- 
ceros ring and securing the metallicity spread in the outer 
thick disk. Moreover, precise predictions would require nu- 
merical simulations with star formation and metal enrichment. 

Finally, we note that the so -called TriAnd clump feature 
(Rocha-Pinto et al. 2003, 2004) is sometimes discussed in as- 
sociation with the Monoceros ring. It is worth emphasiz- 
ing that whether these structures are related remains a matter 
of debate. TriAnd is spatially disjoint from the Monoceros 
ring, with c4eiio > 20 kpc and I = 100- 150° compared to 
<^heiio ~ 12 kpc for the Monoceros stream. Some models sug- 
gest t hat a single disruption event can explain both features 
(e.g.,|Penarrubia et al. 2005). However, it is also possible that 
the TriAnd clump is a disjoint disruption event, or possibly 
a secondary ring feature, of the type colored in green in Fig- 
ure|5] More exhaustive comparisons between data and models 
will be required to definitively settle these issues. 

Following up on the findings reported in Figs. [3] and |4] Fig- 
ure [6] illustrates the effect of each satellite passage S1-S5 on 
the global morphological evolution of the disk. Images were 
created after the disk has relaxed from the encounter with 
each infalling subhalo according to the definition in §2.5. The 
high-surface brightness thin disk component discussed above 
is visible in all edge-on views of the disk. The first satellite 
passage S 1 generates a pronounced warp in the disk beyond 
some radius. The second interaction with the most massive 
subhalo S2 has a dramatic effect on the disk structure causing 
substantial flaring above the disk plane. Subsequent satellite 
accretion (S3-S5) does not substantially perturb the disk fur- 
ther. The disk flare is visible at low surface brightness limits 
below 26 - 27 mag arcsec" 2 and becomes particularly promi- 

8 Note that the inner thick disk may have a different metallicity, so the com- 
parison should be performed with the outer thick disk, which likely formed 
during the same heating events that produced the ring. Note also that if a 
similar ring-like feature is generated in the gas distribution, associated star 
formation and enrichment in the ring itself can enhance the metallicities of 
the stars in the stream. 



nent after the second impact. This fact suggests that only a 
few interactions with satellites with ACDM motivated inter- 
nal properties and orbital parameters can generate the bulk of 
morphological evolution of a stellar disk. We return to this 
point below. 

3.3. Disk Flaring 

Among the most striking morphological features in the fi- 
nal edge-on disk is the pronounced flare, which is particularly 
visible at low surface brightness levels (Figures [3] and |SJ. A 
more quantitative measure of the flare induced by the accre- 
tion events in our simulations is given in Figure [7] This figure 
shows scale height profiles of the disk, z(R), viewed edge-on. 
Profiles are normalized to the initial disk scale height, Zd, and 
are plotted as a function of projected radius in units of the disk 
radial scale length, Rj. 

The flaring of the disk can be formally described by the 
increase of the vertical scale height. For the purpose of our 
analysis, we define the scale height at any given radius to be 
the vertical distance above the disk midplane where the edge- 
on surface brightness falls off from its maximum along the z = 
disk plane by A/i = 1 mag arcsec -2 , /z(z) = /j.q+1 where /iq = 
fJi(z = 0). Though disk scale heights must be formally derived 
by means of fitting the particle distribution to an appropriate 
functional form (e.g., exponential or sech 2 law), scale heights 
as quantified here are well-defined and require no assumptions 
for their interpretation. 

The left panel of Figure|7]shows scale height profiles for the 
initial and final disk. Results are presented using both the def- 
inition of scale height introduced above, /i(z) = {1q+ 1, and that 
for which scale heights are defined to be the vertical distance 
above z = where the edge-on surface brightness falls off from 
the central maximum value by 2 mag arcsec" 2 , /i(z) = /j,q + 2. 
The final disk shows considerable flaring beyond ~ 2 - 3Rd 
indicating that encounters with CDM substructure can be re- 
sponsible for producing notable flares in stellar disks. Re- 
markably, the scale height of the disk near the solar radius has 
increased in excess of a factor of 2 as a result of the interac- 
tions with satellites S1-S6. We checked that the scale height 
of the same disk galaxy evolved in isolation for a timescale 
equal to that corresponding to the final disk (~ 8 Gyr) has 
grown by only about 10% indicating the quality of the initial 
conditions and adequate resolution of the simulations. 

The right panel of Figure|7]shows the evolution of the scale 
height profile of the disk caused by individual satellite pas- 
sages. After the encounter with satellite SI, the inner disk 
(R < 2Rd) appears unaffected by the accretion event despite 
that fact that SI is quite massive, with an initial mass of 
~ 30% of the disk mass, and is characterized by a small peri- 
center (r per i ~ 1 -2Rd). In contrast, the scale height of the outer 
disk (R > 3.5Rd) raises considerably giving rise to a distinct 
flare. The second event, which involves the most massive 
satellite S2 (M su b ~ 0.6Md; s k), generates a stronger flare by 
greatly increasing the scale height of the outer disk. Though 
the scale height has now grown throughout the entire disk, 
the inner disk regions still exhibit substantial resilience to the 
encounter. Indeed, at R = Rj the scale height only raises by 
~ 50% compared to approximately a factor of 4 increase at 
R = 4R d . 

This difference in the robustness of inner and outer disks 
is attributed to the larger binding energy of the former and to 
the presence of a massive, central bulge that acts as a sink of 
the orbital energy of the infalling satellites. The combined ef- 
fect of the remaining satellite impacts (S3-S6) on the structure 
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FIG. 7. — Disk flaring. Scale height profiles, z(R), of the disk viewed edge-on as a function of projected radius in units of the disk radial scale length, Rj. The 
scale height is defined to be the vertical distance from the disk plane where the surface brightness drops by 1 mag arcsec - "" from its maximum along the z = disk 
plane, /i(z) = no + 1 where fig = H(z = 0). Surface brightnesses are averaged about the z = and R = planes and all profiles are normalized to the initial disk 
scale height, z.d- Left: Scale height profiles for the initial (thin lines) and final (thick lines) disk. Solid lines show results obtained using the above definition of 
scale height, while dotted lines correspond to scale heights defined to be the vertical position above z = where the edge-on surface brightness falls off from the 
central value by 2 mag arcsec~ 2 . The initial disk is constructed with a constant scale height explaining why the corresponding curves are flat. The vertical dotted 
line indicates the location of the solar radius, Rq. A conspicuous flare is evident in the final disk beyond ~ 2 — 3Rj as a result of the gravitational interaction 
with CDM substructure. Right: Evolution of the scale height profile of the disk. Various lines correspond to different satellite passages from S 1 to S6 and scale 
heights are measured using the standard definition, /i(z) = fJ-o + 1 ■ The encounters with the first two satellites give rise to a distinct flare in the outskirts of the 
disk. The combined effect of the remaining satellite passages (S3-S6) is much less dramatic indicating that the second, most massive encounter is responsible for 
setting the global disk structure. The inner disk appears much less susceptible to damage by the infalling subhalos owing to its large binding energy and presence 
of sinks of the orbital energy of satellites, such as the bulge component. 



of the disk is much less dramatic increasing the scale height 
only slightly compared to passage S2. This finding suggests 
that subsequent accretion events by already thickened disks 
induce much smaller changes in the disk scale height com- 
pared to th e initial encounters . This effect was also noted pre- 
viously by Oui nn et ail 0993) using numerical simulations in 
which halos were treated as rigid potentials. In broad terms, 
the dynamical and morphological evolution of a galactic disk 
subject to a cosmologically-motivated subhalo merger history 
are driven by the most massive accretion event. 

Given the fact that the self-gravity of the disk grows weaker 
as a function of distance from the center, it is not unexpected 
that the scale height of the final disk should increase with ra- 
dius. In the thin disk approximation, we may consider the to- 
tal vertical energy per unit area of the disk, e z = t z + Wdd + Wds, 
where t z is the disk vertical kinetic energy, and Wdd and 
Wds are the disk potential energy densities associated with 
the disk self-gravity and disk-spheroid gra vity, respectively 
dToth & Ostrikerl 1 19921: iBenson et al"1l2004l) . In this formu- 
lation it is implicitly assumed that the spheroid includes the 
dark matter halo and any bulge component. For a thin disk of 
varying scale height, z(R), we expect 

w dd = GaZ 2 d (R)z(R) (4) 

w ds = G[3^ d (R)Ps(R)z 2 {R). (5) 
Here, p s (R) is the spheroid mass density averaged over the 
disk scale height at radius R and ~E d (R) is the projected disk 



surface density. The constants a and /3 are geometrical factors 
of order unity that will depend on the profile shapes of the disk 
and spheroid components, respectively. 

Now consider the impact of a dark subhalo with the disk. 
Some fraction of the orbital energy of the satellite (~ M^V^, 
where V sat is the orbital velocity), will be delivered as vertical 
kinetic energy to the disk. The disk will v irialize on a dy- 
namic al timescale such that t z ~ Q.5w dd +w d s dTofh & Ostrikerl 
119921) . When a subhalo passes through the disk, it fre- 
quently triggers global modes (e.g., spiral arms, warps, bars) 
which can be very ef ficient in redistributin g its orbital energy 
throughout the disk (Sellwoo d et alj|1998l) . Moreover, if the 
infalling satellites are extended, like in the case of our sim- 
ulations, most disk stars will also directly receive kinetic en- 
ergy during the interaction. Both facts suggest that the orbital 
energy of an accreting subhalo will not be deposited locally 
at the point of impact as it was assumed by Toth & Ostrikerl 
( 1992), but rather globally across the entire disk. 

The simplest assumption for energy deposition is that it is 
roughly constant in radius such that e z (R) — » e z (R) + Ae z , with 
Ae z nearly independent of R. Moreover, due to the fact that 
rotational energy dominates over random motions in the plane 
of the disk, the predominant heating will be in the vertical 
direction. In this limit we can assume that the projected ra- 
dial profiles of the disk and spheroid remain unchanged. The 
global change in vertical energy will demand that the scale 
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height changes as: 
Az(R) oc 



X d (R)[aX d (R) + 2z(R)f3p s (R)Y 



(6) 



In the limit where the disk surface density dominates the 
restoring force, the disk scale height will increase as Az(R) oc 
Sj 2 (/?). The presence of additional sinks of the orbital energy 
of satellites such as a massive, central bulge or a concentrated 
dark matter halo will further act to suppress the increase of 
scale height at small R. Indeed, as demonstrated in Paper II, 
a dense bulge is efficient in reducing the damaging effects of 
infalling subhalos on stellar disks. 

Note also that once the disk is thickened to drive z(R) to 
be quite large, we expect the second term in the denominator 
of eq. |6) to prevent any further vertical puffing of the disk. 
Generally speaking, this suggests that the most massive of the 
past encounters will set the scale height of the old thick disk, 
which is directly confirmed by the findings presented in the 
right panel of Figure [7] 

It is particularly interesting to investigate how the pro- 
nounced flare we predict in our simulations would manifest 
in the surface brightness distribution of an edge-on disk. The 
relevant analysis is presented in Figure [8] This figure shows 
an edge-on "in-plane" (z = 0) surface brightness profile for 
both initial and final disk as a function of projected radius, 
R. The calculation assumes again a stellar mass-to-light ratio 
of M*/L = 3. A steepening of the exponential slope can be 
seen in the profile of the final disk beyond a projected radius 
of R ~ 10 kpc. The entire surface brightness distribution of 
the final disk also decreases as a result of the encounters with 
the accreting satellites. Specifically, the surface brightness in 
the outskirts (R ~ 6Rd) drops from /i ss 27 to \i » 29 mag 
arcsec" 2 . 

Furthermore, the central surface brightness decreases from 
the initial value of fj, ~ 21.2 mag arcsec" 2 by A/i « 0.6 mag 
arcsec -2 . This is consistent with the decrease in the peak of 
the surface brightness profiles at R = between initial and 
final disk seen in Figured] As we discussed earlier, this drop 
as well as the overall decrease of the brightness distribution 
can be explained by the fact that infalling satellites heat stars 
efficiently above the disk plane. 

A low-surface brightness phenomenon seen in many 
disk galaxies is the tendency for exponential d i sks to 
be truncated at faint magnitudes (|van derKruitl fl979l 



pe truncated at tain t magnitudes ( van der ivruit 
Ivan der Kruit & Searld 1 198 it Ide Jong et all l2007bl) . Such 
breaks were originally discovered in edge-on systems and ap- 
pear to be less readily observed in face-on galaxies because 
of the line-of-sight projection and the associated lower sur- 
face brightness. It is important to emphasize that even though 
the behavior of the final edge-on disk depicted in Figure[8]re- 
sembles observed disk truncations, the steepening of the ex- 
ponential slope disappears when the disk is viewed face-on. 
This suggests that the 3D stellar distribution does not actu- 
ally truncate. Thus, the steepening of the surface brightness 
profile p(z = 0,R) reported here does not constitute a separate 
physical effect. On the contrary, it is a consequence of the 
flare. 

In the outskirts, the binding energy of the disk is low and 
stars are efficiently heated above the disk plane by the in- 
falling satellites. As a result, at large distances from the cen- 
ter, very little material remains in the disk plane which gives 
rise to a distinct drop-off in the surface brightness beyond 
R ~ 3.5Rd- On the other hand, the larger self-gravity of the 
inner, exponential disk and the presence of a massive bulge 
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FIG. 8. — Edge-on surface brightness profiles of the stellar disk as a func- 
tion of projected radius, R. The final disk (dotted line) exhibits a steepening 
in surface brightness distribution beyond R ~ 3.5Rj compared to the initial 
disk (solid line). The peak of the surface brightness profile of the final disk 
decreases by A/i r; 0.6 mag arcsec -2 from the initial value. At large dis- 
tances from the disk center a differential decrease in the surface brightness of 
Afi 2 mag arcsec - ' is established. Disk flaring due to infalling CDM sub- 
structures is responsible for generating breaks in the edge-on surface bright- 
ness profiles of galactic disks. 



that acts as a sink of the satellites' orbital energy increase the 
robustness of the inner disk to accretion events and prohibit 
the development of such a break in the exponential slope at 
smaller radii. 

4. DISCUSSION 

In what follows we discuss the implications of our findings 
and caveats of the present study. 

4.1. Implications 

We have demonstrated that encounters with CDM substruc- 
ture generate a wealth of morphological signatures in disk 
galaxies. These include conspicuous flares, bars, low-surface 
brightness ring-like features in the outskirts of the disk, faint 
filamentary structures above the disk plane, and a complex 
vertical morphology that resembles the commonly adopted 
thin-thick disk pr ofile s used in the analysis of disk galaxies. 
iFont et ail d2001l) and iGauthier et al.l d2006) performed sim- 
ilar numerical investigations of the gravitational interaction 
between multi-component disk galaxies and infalling subha- 
los. Both of these investigations considered the population of 
surviving substructures at z = and both studies found negli- 
gible effect on the global structure of the disk. 

In particular, the aforementioned studies showed only mild 
evolution of the disk scale height as a function of radius due 
to the satellite population. In contrast, the present study re- 
ports substantial disk flaring stemming from the gravitational 
interaction with infalling CDM subhalos. As Figure [7] indi- 
cates, the flaring we find is associated with an increase of the 
disk scale height by more than a factor of 2 near the solar ra- 
dius. The primary reason for this discrepancy is that we have 
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followed the formation history of the host halo since z ~ 1, 
and consequently we have considered impacts of satellites 
onto the disk that are significantly more massive tha t those 
incl uded in the numerica l experiments of iFont et al.l d2001l) 
and lGauthier et all (120061) . 

This difference is critical because massive subhalos on 
highly eccentric orbits at early epochs become preferentially 
disrupted or suffer substantial mass loss during their orbital 
evolution, and so the y are more likely to be absent from the z 
I i population (e.g.. | Zentner & Bullockl 120031 : iRravtsov et alj 
l2004tlGao et al.ll2004l) . To capture these relatively short-lived 
objects, it is necessary to trace the interaction history of the 
central halo with nearby, massive systems as we have done 
here. Interestingly, orbital evolution models of the LMC and 
SMC using 3D velocities constrained by recent proper mo- 
tion measurements suggest that these massive Galactic com- 
panions may curre ntly be on their first passage about the MW 
dBesla et al. 2007). If confirmed, this is in line with the idea 
that massive objects would be present today only if they were 
recently-accreted. 

As discussed in the introduction, it is tempting to associate 
the faint structures resulting from the satellite-disk interac- 
tions (Figures [3] and [6) with the recently-discovered complex 
disk phenomena seen in and around the MW and M31, and 
in other nearby and distant dis k galaxies. Stellar counts in the 
MW f rom 2MASS and SDSS (iMomany et al.l2006l:lBlireTaII 
2007) identify stellar streams, many of which may be associ- 
ated with recent merger or fly-by activity from Local Group 
galaxies, or other perturbations in the galactic disk. These 
studies also fin d a great deal m ore structure at radii beyond 
Rqc ~ lOkpc dBell et al.ll2007l) . as well as a many-fold in- 
crease in th e scale height of the d isk, compared to that at the 
solar circle (Momany et al. 20061). 

Figure [6] also illustrates that encounters with infalling 
subhalos can generate warps in the outer disk. Warp exci- 
tation by tidal interactions with small companion galaxies 
constitutes an attractive mechanism. Indeed, it has recently 
been argued that the origin of the Galactic warp can be 
attributed to the interact ion of the MW disk with the LMC 
(Wein berg &Blrt3l2006h . In more distant galaxies, extremely 
deep imaging has long revealed that multiple-component, 
thick-disk, and more extended diffuse structures are not 
limited to MW and M3 1 extending the applicability of our 
results beyond the Local Gro up (e . g.. ISackett et alJ Il994t 
Leaueux et alj|1998t IXbe et alj|1999t iDalcanton & Bernstein! 
2002t IZibetti & Fergusonl J2004t IDalcanton et all 120051: 
Yoachim & Dalcantonll2006l: fde Jong et al.ll2007ai) ~ 

The ring-like features shown in Figure|3]are particularly in- 
triguing in light of the Monoceros ring structure of the MW. 
In Figure [5] we presented a comparison between the ring- 
like structures in our simulations and the Monoceros ring and 
demonstrated that there is a striking agreement in a num- 
ber of properties including the R and z spatial distributions, 
tight velocity widths and line of sight velocity dispersions. 
To this end, the present study suggests that dynamically cold 
rings of this kind are generated naturally as a result of im- 
pacts with halo substructure, which can excite perturbations in 
disk stars. This phenomenon is different from the well-known 
"shell" structures that arise in stars that are accreted from dy- 
namically cold systems that interact with larg er galaxies (e.g., 
iHernquist & Ouinall988fclHelmi et all l2003). In our simula- 
tions, the ring is caused by a perturbation within an initially 
cold disk. 

It is important to emphasize that our simulation set was 



not designed to produce a ring-like structure, nor any of the 
other morphological features observed in the remnant stellar 
disk. Our campaign was designed to mimic the general pro- 
cess of satellite infall in the prevailing cosmology and as such, 
we expect our results to be general consequences of this pro- 
cess as well. In comparison to models that attempt to explain 
the Monoceros stream and extended M31 disk structure us- 
ing minor merger events that directly deposit stars, our sce- 
nario requires no special orbital configuration for the infalling 
satellites. In order to explain the origin of the aforementioned 
features, accreted-star models seem to require fai rly circular, 
prograde, or (at least) low-inc lination orbits (e.g. jHelmi et alJ 
l2003t Penarrub ia et alJ 120051) . In studies of dissipationless 
simulations, orbits of this type are quite rare for subhalos that 
penetrate deep into the host halo and interact with disks (e.g ., 
i Ghiena et all[l998t iKnebe et alJl2004t IZentner et all2005albl : 
lBensonir2005l) . The subhalo orbital eccentricity distribution 
presented in the right panel of Figure Q] confirms this conclu- 
sion. 

More circular orbits may be motivated by numerical exper- 
iments of satellite decay in disk-dominated or flattened sys- 
tems that have repor ted ra pid circularization of th e orbit (e.g., 
l Ouinn & Goodmanlll986t iPenarrubia et alJl2002h iMeza et al] 
2005). However, because dynamical friction depends sen- 
sitively on the mass ratio of the interacting systems, among 
other things, the mass of the satellite must be a non-negligible 
fraction of the disk mass for such mechanism to be viable, and 
in the very least subhalo orbits will not be driven to nearly 
circular prior to an encounter with the disk. The extended 
disk components and morphological features we report in the 
present study arise naturally in the context of a typical ACDM 
accretion history. 

While ring features in the simulated disks capture many 
of the properties of the Monoceros stream, there are a num- 
ber of observational constraints that cannot be directly ad- 
dressed with our simulations and could potentially falsify our 
proposed model. A promising way to distinguish between 
the two different scenarios for the formation of the Mono- 
ceros ring is to establish that the stream contains stellar pop- 
ulations that are not expected to pre-exist in the outer disk. 
Interestingly, there have been indications that some globular 
clusters may be associated with the Monoceros stream (e.g., 
ICrane et alJl2003l) . While this association is still a matter of 
debate, the existence of such globular clusters would point to 
an external origin for the Monoceros structure. 

In principle, the metallicity distribution of stars in the 
stream can also be used to constrain the competing mod- 
els. Indeed, in our scenario the metallicity of the stream 
should be comparable to that of the thick disk stars sur- 
rounding the stream. If the outer thick disk is truly as 
metal rich as the canonical thick-disk value , [Fe/H] w -0.7 
to -1 dGilmore et al.lll995l: llvezic et al.)l2 008). the higher ring 
metall icities estimated by iCrane et all d2003l) and llvezic et al.l 
(2008) would be consistent with our scenario. Taking this 
canonical thick-disk value, t he lower ring metallicities re- 
ported by lYannv et"ail (120031) would seem to be problematic 
for our model. Overall, in order to use metallicity as a re- 
liable constraint on our model, we would have to refine ob- 
servational measurements for the Monoceros ring, secure the 
metallicity spread in the outer thick disk, and perform numer- 
ical simulations with star formation and metal enrichment. 

In the context of distinctive observational signatures pro- 
duced in stellar disks by accreting subhalos, an intriguing re- 
sult is reported in Figure [7] This plot demonstrates that the 
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scale height of the simulated disk increases with radius as a 
result of the gravitational interactions with infalling satellites. 
Indeed, disk flaring appears to be a natural result of encoun- 
ters with CDM substructure and one of the most important 
observational consequences of the present study. Many late- 
type galaxies (including the MW and M31) exhibit a great 
deal of non-smooth "structure" in their stellar distributions, 
often with clear flaring at large galacto-centric radii. 

Actual disk-flaring, such as that described in this 
paper, is seen in the MW in both the stellar disk 



p . 

dLopez-Corredoira et al.1 120021; iMomanv et all 2006) 
and atomic hyd r ogen layer (e!g^ iMerrifieldl 11992 : 
iNakanishi & Sofud 120031) . Flaring is also observed in 
external galaxies seen ed g e-on, in their s t ellar l ight (e.g., 
Ide Griis & Pelefieri 119971: iNaravan & Joel |2002|) (th ough 
more-rarely in lo wer-mass galaxies; de Griis & Peletie rll 1 997t 



Seth et alj|2005l). a s well as in Hi gas (e.g., | Brinks & Burtonl 



1984t lOllindl 19961: iMatthews & Wood 2003). Extending the 



search for disk flaring in external systems would provide 
a new observational test for the ACDM paradigm and its 
competitors on non-linear, galactic scales, and thus offer 
unique opportunities to constrain cosmological models. 

The structure and kinematics of galactic disks have long 
been recognized as a key constraint on structure forma- 
tion models. Our results have intriguing implications for 
tests of the nature of dark matter through detailed ob- 
servations of galactic structure. Alternatives to the stan- 
dard CDM paradigm such as mo dels with warm dark mat- 
ter or decaying da rk ma tter fe.g.. iHogan & D alcanto dl2000l 



Avila- Reese et al.l 1200 It IC embrano s et al.l 12005b Kaplinghat 
20051; [Strigari et al. 200 7b|). or non-st andard power spec- 
tra (e.g.. lKamionkowski & Lid dle 2000) predict substantially 
different accretion histories for galaxy-sized dark matter halos 
as w ell as different satellite ab undances and structural proper- 
ties dZentner & Bu llock 200j). 

We have assessed the sensitivity of disk flaring to the den- 
sity structure of the infalling systems by adopting satellite 
models with radically different density profiles, but requir- 
ing that the total bound mass of each object be exactly the 
same. Specifically, we repeated the first two satellite passages 
SI and S2 modeling the infalling subhalos using a constant 
density core (p(r) oc r" 7 , where 7 = 0) and adopting the same 
outer, /?, and intermediate, a, slopes of their CDM coun- 
terparts. The size of the density core was chosen to be ap- 
proximately equal to 10% of the tidal radius of each subhalo 
(> 2 kpc). 

We find that cored satellites of equal bound mass produce 
much smaller amounts of disk flaring at the solar radius and 
at larger radii, and smaller overall disk morphological evolu- 
tion compared to the standard CDM cosmological satellites 
whose density distribution is well described by a cusp slope 
of p(r) oc r~ l . This is not unexpected since the former are 
more susceptible to mass loss and will, on average, approach 
and penetrate the disk with smaller amounts of orbital energy 
available to be converted into random motions of disk stars. 
These findings suggest that differences in the properties of the 
subhalo populations should be imprinted on the fine structure 
of disk galaxies. 

Such differences can be used to provide fundamental con- 
straints on structure formation theories by distinguishing be- 
tween competing cosmological models. Indeed, specifically- 
designed surveys such as SEGUE and RAVE and planned rev- 
olutionary instruments like SIM and GAIA will be able to 
quantify the structure and kinematic properties of the MW 



disk. Moreover, future surveys like SEGUE-2, PanSTARRS, 
and LSST will provide ever more precise star-by-star maps 
of the Galaxy and multi-object spectrographs on thirty-meter 
class telescopes will enable star count maps for many galaxies 
in the local volume, similar to the maps being made for M3 1 
and M33. Confronting the resulting data sets with theoretical 
predictions such as those of the present study will allow to test 
whether the detailed structure of galactic disks is as perturbed 
as predicted by the CDM paradigm. 

We expect disk morphological evolution in galaxy-sized 
dark matter halos, M vlr ~ 1O 12 M , to be driven by the few 
most massive satellite accretion events M su b ~ 10 10_11 Mq ~ 
A^disk- Indeed, mergers of this size dominate the mass buildup 
in galaxy-sized CD M halos (e.g., Zentner & Bullocki 120031 ; 
iPurcell et al.l 120071) . In addition, massive accretion events of 
this kind should have been relatively common since z ~ 1 . We 
have showed this directly using our four well-resolved cosmo- 
logical jV-body halo s in § 2.2. This result is corroborated by 
Stewa rtet al.l J2007) who used thousands of halos extracted 
from a ACDM Af-body simulation (erg = 0.9) to show that 
~ 80% (~ 70%) of MW-sized objects should have accreted 
at least one subhalo with mass comparable to that of the disk 
in our controlled simulations, M su b > 0.02M v ; r (> 0.05M v j r ) 
since z ~ 1 . Moreover, ~ 60% should have accreted a sys- 
tem with more than twice that mass over the same period 
(> 0.1M v i r ). Of course, the importance of these accretion 
events for both driving and altering disk morphologies will 
depend on the orbital evolution of the satellites once accreted 
into the host halo virial radius. However, given that dynami- 
cal friction will drive massive subhalos towards the central re- 
gions of the host halo within a few Gyr, these results serve as 
general motivation that significant central encounters should 
be fairly commonplace in the history of galaxy-sized halos. 

The above discussion is particularly relevant in the context 
of the survival of thin, stellar disks to hierarchical satellite ac- 
cretion. Figures[3] |4]and |6]illustrate that while infalling satel- 
lites substantially perturb stellar disks affecting their overall 
morphology, a significant thin disk component can survive 
bombardment by CDM substructure that is an appreciable 
fraction of the disk mass. Though most MW-sized galaxies 
are expected to have accreted an object at least as massive as 
the disk itself since z ~ 1 (Figure[T]i we have explicitly ignored 
these violent interactions in the present study in an attempt to 
be conservative. This is because our main goal was to exam- 
ine the formation of distinct morphological features about a 
well-preserved disk galaxy similar to the MW and M3 1 . An 
important future study will be to model the impact of a satel- 
lite as massive as (or more massive than) the disk itself in or- 
der to test the robustness of disks in this general framework. 
We intend to examine this issue fully in a forthcoming work. 
The results presented here provide some encouragement that 
such an experiment may not prove ruinous to the survival of 
thin disks in a ACDM universe. 

Of specific concern for hierarchical models of structure for- 
mation is t he prevalenc e of old (~ 10 Gyr) s tars in the thin 
disk (e.g., IWvsd |2001|) . lOuillen & Garnettl d2000l) showed 
that the age-velocity dispersion relation in the solar neigh- 
borhood for disk stars was relatively flat between ~ 2 and 
~ 10 Gyr, with a distinct jump in the velocity dispersion for 
older stars, as perhaps would be expected if thick disk forma- 
tion were driven by a single an cient event. However, u tilizing 
a sample ~ 75 times as large, Nordstr om et alJ (|2004) found 
a continuous increase of stellar velocity dispersion with age, 
and argued that ongoing heating could provide an explana- 
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tion. ISeabroke & Gilmorel (120071) used the same data set and 
similarly concluded that vertical disk heating does not satu- 
rate at early times. Regardless, the age distribution of stars 
in the MW suggests that a significant fraction of the thin disk 
(a z < 20 kms" 1 ) was in place by z ~ 1. 

It is tempting to argue that the existence of an old, thin stel- 
lar disk would imply an absence of satellite impacts over a 
period that corresponds to at least the age of the oldest stars in 
the thin disk. If an impact occurs and destroys the disk, then 
a thin disk component must be generated later by secondary 
processes such as gas inflow. This would suggest that thick 
disks originating during damaging encounters with substruc- 
tures must be invariably older than the oldest stars in the thin 
disk, making thin and thick disks distinct in age as well as 
scale height. While this is true for the MW and other external 
galaxies, our results are consistent with a picture in which a 
thick disk arises naturally in the context of a typical ACDM 
accretion history while at the same time a substantial thin disk 
component survives the violent gravitational encounters with 
CDM substructure. Our results thus indicate that the absence 
of thick disks in significant fraction of disk galaxies, rather 
than ubiquitous thick disks, would represent a challenge for 
the CDM model. 

The findings of the present study have implications for the 
formation of thick disks whic h constitute a ubiquitous com- 
ponent of disk galaxies (e.g ., iDalcanton & Bernstein! 120021 : 
lYoachim & Dalc anton 2006). Several formation mechanisms 
have been proposed to explain their origin and properties 
with a growing body of theoretical and observational ev- 
idence favoring a merger-driven scenario. There are two 
qualitatitavely different mechanisms that are associated with 
this formation scenario. In the first, thick disk stars form 
initially in the thin disk and then are dynamically heated 
to lar ge scale heights by encounters with orbiting satellites 
(e.gJOuinn et alJl993l:IWalker et alJl996HRobin et all 19961: 
I Velazquez & Whitelll999t IChen et all 1200 lh . In the second, 
thick disk stars form in external galaxies and subsequently 
are deposited by accretion events at large scale heights (e.g . , 
IStatlerll 1 988t I Abadi et alll200H lYoachim & Dalcantonl2005l) . 

The results presented in this paper suggest that at least 
part of a galaxy's thick disk component may plausibly orig- 
inate from the gravitational interaction between an existing 
thin disk and infalling satellites with mass functions, density 
structures, and orbital distributions of the kind expected in the 
ACDM paradigm of structure formation. CDM substructure 
increases the scale height of stellar disks (Figures [4] and 17) 
and should be regarded as being at least partially responsible 
for the origin of thick disks. Indeed, this conclusion is sup- 
ported by observ ational studies of the vertical distribution of 
stars in the MW (IChen et all 1200 lh and star count dat a from 
a number of Galactic sample fields (Ro bin et al.l fT996). Pho- 
tometric constraints in combination with kinematic data from 
these investigations favor the mechanism in which the thick 
disk of the MW formed through the heating of a preexisting 
thin disk, with the heating mechanism being the merging of 
a satellite galaxy. Additional supporting proof comes from 
recent HST studies of resolved stellar populations in nearby 
edge-on disks fin ding evidence for continuous vertical heating 
dSeth et al.l2 005) and chemical abundance analysi s of F and G 
dwarf stars in the Galactic thin and thick disks dBensbv et alj 
12001 . 

In contrast to the previous investigations, other "thin-thick" 
disk kinematic studies favor models in which the thick disk 
forms from direct accretion of stars from infalling satellites. 



For example. lYoachim & Dalcanton] d2005l) presented GMOS 
kinematic measurements of edge-on galaxies FGC 1415 and 
FGC 227 and found that both thick disks were rotating much 
slower compared to the thin ones. In the case of FGC 227, 
the thick disk even showed weak signs of counterrotation. As 
argued by these authors, the detection of very slowly rotat- 
ing or counterrotating thick disks support an accretion origin 
for thick disk stars. While it is sufficiently unlikely that the 
thin disk forms from subsequent accretion of gas with angular 
momentum opposite to that of the original disk, disks heated 
by satellites do exhibit a v ertical gradient in rotational veloc- 
ity (lHayashi & C hiba 2006) which is co nsistent with what is 
inferr ed for the thick disk of the MW (lAllende Prieto et all 
2006). Interestingly, analysis of the mean azimuthal velocity 
of disk stars at the solar radius in our simulations also reveals a 
vertical gradient in rotational velocity of ~ -20 kms" 1 kpc" 1 
between 1 and 3 kpc from the disk plane, which is in good 
agreement with the result s of dHavashi & C hiba 2006|) and 
dAllende Prieto et alJl2006l) . 

In addition to slower rotation, MW thick disk stars ex- 
hibit larger vel ocity dispersion s compared to stars in the 
thin disk (e.g., IChib a & Beers] l2000h . This is relevant as 
our simulations do reveal that the stellar disk becomes 
substantially heated in all three directions by the encoun- 
ters with CDM substructure. Indeed, as we show in Pa- 
per II, the disk velocity ellipsoid at the solar radius, R Q = 
8 kpc, increases from (<7s,<70,<7 z ) = (31,24, 17) kms -1 to 
(aR,a<p,a z ) ~ (61, 49, 31) kms" 1 . For reference, the veloc- 
ity ellipsoid of the thick di sk of the MW is estimated to 
be - (46, 50, 35) km s" 1 by [Chiba & Beers! d2000h and ~ 
(63, 39, 39) kms" 1 by ISoubiran et alJ d2003l) . However, it is 
important not to over-interpret this comparison because the 
final dispersions are calculated considering all stars instead 
of using only those that belong to the thick disk component 
according to the decomposition analysis of § 3.1. 

Chemically, thick disk stars in the MW are more metal- 
poor than stars in the thin disk (e.g., Chi ba & Beers! l2000l) . 
are significantly enhanced in a-elements compared to thin 
disk stars of similar iron abundanc es which suggests a n ex- 
tended star formation history (e.g., iBensbv et alj|2005|). and 
exhibit no vertical metall icity gradient ( Bensby et al J 120051 : 
lAllende Prieto et all 120061) . Most of the above properties 
have been also confirmed in external systems (e.g. JSeth et alj 
120051) . Obviously, our dissipationless simulations can neither 
verify nor disprove any of the aforementioned trends. How- 
ev er, relevant in this conte xt are the collisionless simulations 
of lHavashi & Chibal d2006l) reporting that satellite-disk inter- 
actions constitute an efficient mixing mechanism capable of 
erasing any metallity gradient present in the original thin disk. 

In summary, there is no definitive observational or theoret- 
ical evidence to conclusively rule out either of the merger- 
driven scenarios for the origin of thick disks. To this end, 
the existence of very slowly rotating or counterrotating thick 
disks in a significant fraction of disk galaxies would be prob- 
lematic for the "vertical heating" model. Most likely, both 
mechanisms do operate at a different degree in forming the 
thick di sk of a galaxy. Thi s is further corroborated by the re- 
sults of lAbadi etaT] (120031) reporting that only ~ 60% of the 
thick disk in their galaxy consisted of the tidal debris of dis- 
rupted systems. Our study demonstrates that thick disks with 
realistic properties can result from the dynamical heating of 
pre-existing thin disks by CDM substructure. More detailed 
modeling of the thick disk properties would require the inclu- 
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sion of gasdynamics, star formation, and metal enrichment. 
4.2. Caveats 

Certainly our approach is not without drawbacks. The most 
evident one is that the present study used only four dissipa- 
tionless galaxy-sized dark matter halos in order to derive "typ- 
ical" orbits and merger histories, and the controlled simula- 
tions of satellite-disk interactions were based on the accretion 
history of just one of these host halos. It will be important to 
explore a range of accretion histories and orbital distributions 
of infalling objects using a larger sample of halos from cos- 
mological A^-body simulations. Moreover, since the orbital 
evolution of substructure once accreted into host halos will 
be influenced by baryonic processes, a self-consistent, hydro- 
dynamical simulation of disk galaxy formation would be re- 
quired to fully refine the conclusions presented here. Unfor- 
tunately, hydrodynamical simulations that produce statistical 
samples of realistic disk galaxies at z = do not yet exist. 

Despite the aforementioned limitations, several facts sug- 
gest that our study has captured to at least a certain extent 
the amount of global morphological transformation that CDM 
satellites can cause to stellar disks. First, all four of the host 
halos we explored showed similar numbers of substantial cen- 
tral mergers, and their accreti on histories are typi cal of sys- 
tems in this mass range (e.g., We chsler et al.ll2002h . Second, 
we used a conservative subset of one of the merger histo- 
ries to seed the controlled satellite-disk encounter simulations, 
and neglected interactions with extremely massive subhalos 
that could prove ruinous to disk survival. If anything, the 
morphological response of most disks should be more pro- 
nounced compared to that we find here. Third, our numerical 
experiments took into account the effects of a realistic bary- 
onic component on the orbital evolution of the infalling satel- 
lites by including a live central galaxy. Finally, most of the 
observational signatures on stellar disks (Figures [3] H] [6] [7] 
and [8) were generated by the single most massive encounter 
(Mmb ~ 0.6Mdi s k, Aid ~ 20 kpc, and r per i < 10 kpc). We iden- 
tify accretion events of this kind in all four of the host halo 
histories studied and, as discussed above, it is reasonable to 
expect that most MW-sized halos should have accreted nu- 
merous such systems since z ~ 1. 

A second caveat is that we have neglected stellar compo- 
nents in the infalling satellites modeling the latter as pure dark 
matter subhalos. Thus, the present study is inadequate to ad- 
dress the contribution of satellite stars to the formation of the 
thick disk and elucidate their ability in driving existing or cre- 
ating their own morphological signatures. Notwithstanding 
this limitation, the present work does demonstrate that dark 
matter subhalos can be capable of generating a wealth of mor- 
phological features in the simulated disks that resemble those 
seen in real galaxies. A third drawback is related to the fact 
that our models track only stars that were in place in a disk at 
Z ~ 1, and therefore apply most directly to old disk stars. We 
have also adopted a primary disk galaxy model motivated by 
the present-day structure of the MW. A more complete inves- 
tigation would have to include the ongoing formation of the 
disk. 

Furthermore, we have only addressed the gravitational in- 
teraction of disk galaxies and dark matter substructures in 
the collisionless regime, a choice with obvious limitations. 
Spiral galaxies contain atomic and molecular gas which can 
absorb and subsequently radiate away part of the orbital en- 
ergy deposited by the sinking satellites. Owing to this prop- 
erty, the heated gas component will eventually resettle to form 



a new thin disk consisting of a younger stellar population 
compared to that of the thick disk. Modeling the evolution 
of a dissipative component in the primary galaxy would be 
required to determine the extent to which the presence of 
gas can alter the dynamical effects of substructure on stel- 
lar disks. Given that the gas fraction in massive disk galax- 
ies at z ~ is typically low (~ 10% of the stellar disk mass; 
iRead & Tre ntham 2005), the role of a dissipative component 
in stabilizing the disks against the violent gravitational en- 
counters with satellites may not be important, but this will 
depend sensitively on whether the gas content was higher in 
the past. Indeed, high gas-fractions at early times may be 
required to expl ain the formation of di sk galaxies in a hierar- 
chical universe (Robertson et al. 20061). 

Finally, the present-day structure of galactic disks origi- 
nates from a complex interplay of effects and a full expla- 
nation requires detailed knowledge of their star formation 
history and chemical evolution, amongst others. In fact, 
the thin and thick disk stars of the MW represent distinct 
populations not only kinematically but also chemically (e.g., 
IWvse & Gilmordll995llBensbv et al.ll2005l) . Adding star for- 
mation as a further ingredient will provide the opportunity to 
track ongoing star formation in an ever re-forming thin disk 
and refine the conclusions presented here. All of these effects 
will be discussed in forthcoming works. 

5. SUMMARY 

Using high-resolution, fully self-consistent dissipationless 
Af-body simulations we have examined the cumulative effect 
of substructure impacts onto thin disk galaxies in the context 
of the ACDM paradigm of structure formation. The main 
goals of the present study have been to (1) demonstrate that 
interactions with CDM substructure do lead to morphological 
signatures and substantial changes in the structure of galactic 
disks and (2) study generic features of the evolution of disk 
galaxies subject to a ACDM-motivated satellite accretion his- 
tory. 

Our simulation campaign utilizes cosmological simulations 
of the formation of galaxy-sized dark matter halos to derive 
the properties of subhalo populations and controlled numer- 
ical experiments of consecutive satellite encounters with an 
initially-thin disk galaxy. The present study extends and ex- 
pands upon past numerical investigations into the interaction 
between disk galaxies and merging satellites in at least three 
major respects. 

First, we incorporate a model in which the infalling satel- 
lite populations are truly representative of those accreted and 
possibly destroyed in the past, instead of the z = surviving 
substructure present in a CDM halo. Second, the primary disk 
galaxy models are fully self-consistent and flexible enough to 
permit detailed modeling of actual galaxies. Third, individ- 
ual satellites are initialized according to the internal structure 
of subhalos culled from high-resolution cosmological simu- 
lations. All of these elements inherent in the modeling en- 
sure that the numerical experiments of satellite-disk encoun- 
ters were performed with a high degree of realism. 

In contrast to what can be inferred from statistics of 
the present-day substructure populations in galaxy-sized 
CDM halos, encounters between massive subhalos (M su b > 
0. 2Mdi s k) and disks should be common occurrences in stan- 
dard ACDM since z ~ 1. Most of the difference is caused 
by the fact that satellites that pass closest to the halo cen- 
ter at early times are precisely those that are most likely to 
suffer significant mass loss or become tidally disrupted. As 
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subhalos on highly eccentric orbits lose mass, the fraction of 
massive satellites with small orbital pericenters declines with 
redshift so that few remain by z = 0. The specific merger 
history we studied involved 6 accretion events of satellites 
with masses, orbital pericenters, and tidal radii of 7.4 x 10 9 < 
M snh /M Q < 2 x 10 10 , r P eri < 20 kpc, and r tid > 20 kpc, re- 
spectively, over a ~ 8 Gyr period. The morphological im- 
pact of these events on an initially-thin (zd = 0.4 kpc) disk 
galaxy (M^k ~ 3.5 x 10 10 M Q ) included the following signa- 
tures, many of which resemble features seen in the MW and 
M31, and in other disk galaxies. 

• The development of a "thick" disk component and the 
survival of a significant "thin" disk component, as char- 
acterized by standard "thin-thick" disk decomposition 
analysis. 

• The growth of a strong bar. 

• The formation of a pronounced flare at intermediate and 
large radii, particularly visible at low surface brightness 
levels. 

• The production of long-lived, low-surface brightness, 
dynamically cold ring-like features in the outskirts of 
the disk similar to observed coherent stellar structures, 
such as the Monoceros ring in the MW. 

• The generation of faint filamentary structures that de- 
velop above the disk plane and (spuriously) resemble 
tidal streams in configuration space. 

Finally, the results presented in this study of a "typical" halo 
merger history highlight the potential difficulty in explain- 
ing the fraction of observed featureless disk galaxies (e.g., 
those with extremely thin discs, non-flared and/or non-warped 
and/or non-barred disk galaxies etc) in the context of the 
CDM model of structure formation. Our findings also raise 
questions regarding whether the existence of pure exponen- 
tial discs or "antitruncated" disks which exhibit an excess of 
surface brightness at very large radii (e.g. JPohlen et al.ll2.007h 
can be explained in a model where mergers are as common as 
predicted in ACDM. 

A relevant intriguing issue concerns the survivability of bul- 
geless disk galaxies. To this end, we have performed numer- 
ical experiments of interactions between satellites and bulge- 
less disk galaxies and found that the latter experience sub- 
stantially more damage compared to their counterparts with 
bulges (see also I Velazquez & Whitdfl999l for a similar con- 
clusion). In this context, the existence of a large number of 
very thin bulgeless disk galaxies would be difficult to recon- 
cile with ACDM. In terestingly, all systems in the sample of 
bulgeless galaxies of Yoachim & Dalcanton (2006) have pro- 
nounced thick disks and there are no signs of companion sys- 
tems in the vicinity of the prototype thin bulgeless disk galaxy 
M33. Directly relevant to this issue are recent, complemen- 
tary studies of edge-on disk galaxies using the Sloan Digi- 
tal Sky Survey database which revealed a significant fraction 
of "super-thin" bulgeless disks with much larger axial ratios 
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compared to typical disk galaxies (Ka utsch et al.ll2.Q06b . On 
the other hand, most MW-sized halos are expected to have ac- 
creted at least one subhalo as massive as the di sk in our con- 
trolle d simulations since z ~ 1 (see FigurefTland lStewart et alj 
l2007h . 

The present work tackles the focused problem of identi- 
fying generic features in disk galaxies that arise in typical 
ACDM merger histories. To this end, we have demonstrated 
that gravitational interactions between existing thin galactic 
disks and CDM substructure should be at least partially 
responsible for the formation of thick disks, disk flares, 
bars, low-surface brightness ring-like configurations and 
faint filamentary structures around disk galaxies. The 
fact that many of these features appear ubiquitous in real 
galaxies is encouraging for the ACDM paradigm of structure 
formation. For ex ample, bars are present in about 70 % of 
disk galaxies (e.g., lKnapen|[l999t |Eskridge et alj|2000l) and 
the occurence of thick disks is even more frequent (e.g ., 
iDalcanton & Bernstein! 120021: lYoachim & Dalcantonl 2006). 
Detailed comparisons with the observed populations of disk 
galaxies would require an extensive series of numerical 
experiments to fully sample the statistical variation in halo 
accretion histories predicted in ACDM as well as an equally 
careful statistical accounting of the fraction of galactic disks 
in the Universe that are indeed featureless. 
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